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ABSTRACT 


Figure  1,  Semi-autonomous  Path  Planning  and 
Obstacle  Avoidance  Research  Results 
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The  Research  Institute  for  Precision  Guided 
Systems  (RIAPGS)  was  awarded  a  21  Month 
grant  on  1  March  07  to  conduct 
groundbreaking  research  on  Agile  Autonomous 
Munitions  (AAM).  The  research  directly  supports 
the  mission  needs  of  the  Air  Force  Research 
Laboratory  Munitions  Directorate  (AFRL/RW)  in 
developing  technology  for  semi-autonomous 
micro  air  vehicles.  Research  supported 
under  the  FY03  RIAPGS  grant  poised  the 
UF-REEF  and  its  partners  to  assume  the  role 
of  a  world  leader  in  AAM  research.  The 
recently  completed  effort  documented  in  this 
final  report  has  built  upon  the  results 
obtained  from  AFOSR  funding  of  the  2003, 
three-year  RIAPGS  grant  and  its  two 
extensions  (8  month  extension  concluding  30 
Nov  2006  and  3  month  extension  ending  28 
Feb  07),  and  leverages  the  capabilities 
established,  in  large  part,  by  AFOSR  DURIP 

grants  awarded  to  the  University  of  Florida  Research  and  Engineering  Education  Facility  (UF- 
REEF).  A  team  of  experts  from  the  UF-REEF,  UF-Gainesville,  University  of  Alabama,  and 
University  of  Michigan  conducted  the  research.  The  effort  was  conducted  under  the  auspices  of 
the  Agile  Autonomous  Munitions  Center  of  Excellence  (AAM  COE),  a  partnership  between  the 
University  of  Florida  (REEF  and  Gainesville  campuses)  and  AFRL/RW.  The  AAM  COE,  which 
supports  the  UF-REEF  Research  Institute  for  Autonomous  Precision-Guided  Systems,  serves  as 
the  focal  point  for  the  research,  development,  design,  analysis,  testing  and  transition  of 
technologies  necessary  to  attain  agile,  smart  munitions  operating  autonomously  or  cooperatively 
in  complex,  uncertain,  adversarial  environments  (see  Figure  1).  The  recently  completed  (30  Nov 
08)  21  month  grant  had  objectives  and  milestones  in  years  2007  and  2008  under  the  following  5 
tasks: 

•  Task  1:  Full  Vehicle  Virtual  Prototyping:  Control  Synthesis  and  Aero-structural 
Characteristics  of  Micro  Air  Delivered  Munitions  (MADM)  and  Micro  Air  Vehicles 
(MAV) 
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2500  500 


•  Task  2:  Integrated  Visual-Servo  Control  and  Path  Planning  for  Unmanned  Air  Vehicles 

•  Task  3:  Aerodynamic  and  Structural  Characterization  of  Flexible  and  Morphing  MAVs 

•  Task  4:  Control  of  Biologically-Inspired  Morphing  for  Variable-Geometry  MAVs 

•  Task  5:  Computational  Aerodynamics  of  Flexible  and  Flapping  Wing  for  MAVs 
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Task  5  was  a  9  month  tasks  that  ended  on  30  Nov  07.  The  remaining  three  tasks  continued  until 
30  Nov  08. 


REVIEW  TECHNICAL  ACCOMPLISHMENTS  BY  TASK 

Task  1:  Full  Vehicle  Virtual  Prototyping:  Control  Synthesis  and  Aero-structural 
Characteristics  of  Micro  Air  Delivered  Munitions  (MADM)  and  Micro  Air  Vehicles  (MAV) 

MILESTONES/STATUS: 

Months  1-9 

•  Complete  derivation  of  methods  for  characterizing  MAV  nonlinear  dynamic  response 
and  aerodynamic  models  for  MAV  configurations.  (Completed) 

•  Static  mounted  MAV  models  will  be  aerodynamically  characterized  and  their 
aerodynamic  coefficients  will  be  modeled  and  provided  to  the  Visualization  Lab  so 
closed  loop  testing  of  the  visual  controller  technology  can  be  simulated  with  innovative 
MAV  aerodynamic  designs  such  as  morphing.  (Completed  at  UF-REEF  Labs  and  UF- 
Gainesville  Labs) 

•  Inner  loop  stabilization  of  MAV  using  visual-servo  control  techniques  derived  in  Task  2 
via  implementation  and  testing  in  the  hardware  in  the  loop  environment.  (Completed) 

•  Derivation,  development  and  implementation  of  virtual  prototyping  methods  associated 
with  path  planning  in  the  image-space  as  derived  in  Task  2.  (Completed) 

Months  10-21 

•  Hardware-in-the-loop  testing  and  evaluation  of  3-dimensional  visual-servo  control 
derived  in  Task  2.  (Completed) 

•  Hardware  -in-the  -loop  testing  and  evaluation  of  geometry  and  obstacle  estimation-based 
MAV  vehicle  control.  (Partially  Completed-  Simulations  run  on  desk  top  computer  and 
not  in  real  time  on  a  field  deployable  laptop) 

GENERAL 

Under  the  leadership  of  Dr.  Rogacki  and  Dr.  Albertani  there  was  significant  progress  made  in 
delivering  a  hardware-in-the-loop  rapid  prototyping  capability  at  the  UF-REEF.  Embedded 
systems  are  designed  to  control  complex  plants  such  as  land  vehicles,  satellites,  spacecrafts, 
unmanned  aerial  vehicles  (UAVs),  aircrafts,  weapon  systems,  marine  vehicles,  and  jet  engines. 
They  generally  require  a  high  level  of  complexity  within  the  embedded  system  to  manage  the 
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complexity  of  the  plant  under  control.  Hardware-in-the-Loop  (HIL)  simulation  is  a  technique 
that  is  used  increasingly  in  the  development  and  test  of  complex  real-time  embedded  systems. 

The  purpose  of  HIL  simulation  is  to  provide  an  effective  platform  for  developing  and  testing 
real-time  embedded  systems.  HIL  simulation  provides  an  effective  platform  by  adding  the 
complexity  of  the  plant  under  control  to  the  test  platform.  The  complexity  of  the  plant  under 
control  is  included  in  test  and  development  by  adding  a  mathematical  representation  of  all 
related  dynamic  systems.  These  mathematical  representations  are  referred  to  as  the  “plant 
simulation.’’  The  UF-REEF  has  established  a  HIL  simulation  capability  for  MAV’s  that  will 
lead  to  a  rapid  prototyping  capability. 

The  Task  1  research  plan  included  testing  the  visual  based  control  system  running  on  the  cluster 
PCs  in  the  REEF  Visualization  Lab  by  flying  a  relatively  stable  off  the  shelf  slow  flier  (similar  to 
MAV  flight  regime)  in  the  small  REEF  wind  tunnel  and  later  in  the  large  wind  tunnel.  The 
general  layout  of  the  test  bed  is  illustrated  in  Fig.  2. 


Figure  2.  Lay-out  of  the  HIL  experimental  test  bed  showing  the  vision- 
based  autopilot  (grey)  and  the  aircraft  systems  (yellow)  in  the  wind  tunnel. 


AIRPLANE  DESCRIPTION 

The  airplane  chosen  for  the  experiment  was  a  biplane  that  had  a  wingspan  of  68.4cm  and  a  total 
wing  area  including  both  wings  of  2642  cm2.  Both  wings  are  flat  plates  and  made  from  thin 
foam  sheets.  As  can  be  seen  in  Figure  3,  the  structure  was  stiffened  using  carbon  fiber  stiffening 
rods.  The  airplane  is  considered  a  profile  airplane  because  it  does  not  have  a  traditional  fuselage. 
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Instead  it  has  two  flat  foam  pieces  shaped  like  a  fuselage  that  are  glued  together  perpendicular  to 
one  another.  Figure  3  is  a  picture  of  the  airplane. 


Figure  3.  Front  and  side  view  of  airplane. 


The  airplane  was  setup  with  conventional  controls  that  consist  of  a  rudder,  elevator,  and  ailerons. 
A  pitot  tube  was  installed  on  the  right  wing  to  measure  free  stream  velocity.  Because  of  the 
fuselage  design,  a  carrying  tray  was  added  to  the  front  of  the  airplane  and  was  used  to  mount  the 
autopilot.  A  ball  joint  was  used  to  join  the  airplane  to  the  mount. 

Communication  with  the  airplane’s  onboard  autopilot  was  conducted  through  the  virtual  cockpit 
software  that  came  with  the  autopilot.  The  computer  on  which  the  software  was  installed  used  a 
communication  box  to  send  signals  to  the  airplane  autopilot.  The  transmitter  was  used  to  trim 
the  airplane  before  maneuvers  and  also  to  input  manual  commands  if  needed.  The  3-DOF  rig 
design  illustrated  in  Figure  4  was  used  to  install  the  Tensor  4-D  aircraft  in  the  wind  Tunnel. 


Figure  4.  View  of  the  airplane  on  the  three  D.O.F.  joint  in  the  wind  tunnel. 

Flight  data  were  recorded  on  the  ground  station  PC  using  the  onboard  telemetry  system.  The 
ground  station,  located  outside  the  wind  tunnel,  was  connected  with  the  cluster  PCs  in  the  REEF 
Visualization  Lab.  The  ground  station  received  data  from  the  aircraft  on-board  system  and  sent 
the  information,  through  an  internet  connection,  to  the  Visualization  Lab. 
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EXPERIMENTAL  SET-UP 


The  hardware  used  in  this  project  includes  the  remote  controlled  airplane,  a  ground  station 
computer,  a  camera,  and  a  server  cluster  used  to  create  a  virtual  environment.  The  scope  of  this 
project  was  to  create  a  communications  system  for  these  four  components  of  hardware  by 
writing  software,  which  interpreted  data  between  them  since  the  camera  and  the  UAV  as  well  as 
the  ground  station,  and  server  cluster  were  in  separate  labs.  The  end  product  software  developed 
enabled  communicated  between  these  systems.  The  quality  of  the  data  transfer  can  be  verified 
by  plotting  control  inputs  and  aircraft  flight  data. 

COMMUNICATION  SOFTWARE 

Two  software  packages  were  designed  and  used  for  two  computers  in  the  HILS  system.  The  first 
program  (Program  1)  is  loaded  onto  the  ground  station  computer  in  the  Wind  Tunnel  Laboratory 
(WT  Lab)  and  the  second  (Program  2),  which  consists  of  a  simple  test  socket  of  a  server,  is 
loaded  onto  the  Visualization  Laboratory  (Vis  Lab)  computer.  These  programs  implement  the 
following  procedures: 

Program  1 

♦  Connect  to  Autopilot  via  a  socket  connection  as  client  using  Virtual  Cockpit 

♦  Connect  to  Vis  Lab  via  socket  connection  as  client 

♦  Retrieve  data  at  6Hz  from  the  autopilot 

♦  Package  data  in  vector  form  containing  the  heading,  pitch,  and  roll.  These  are  2-byte 
integers  with  units  of  radians*  1000) 

♦  Send  data  at  6Hz  to  Vis  Lab 

♦  Receive  data  at  6Hz  from  Vis  Lab 

♦  Output  the  sent  and  received  data  as  well  as  elevator  and  aileron  angles  to  the  output  file 
on  hard  drive  for  graphical  analysis 

♦  Display  transfers  in  the  GUI  on  the  desktop  (specifically  sent  data,  received  data,  elevator 
angle,  and  aileron  angle) 

Program  2 

♦  Connect  to  WT  Lab  via  a  socket  connection  as  a  server 

♦  Receive  data  from  WT  Lab 

♦  Display  data  to  command  prompt 
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RECORDED  FLIGHT  DATA 


Pitch,  roll,  and  yaw  angle  data,  pitch,  yaw,  and  roll  rates,  elevator,  aileron,  and  rudder  servo 
displacement  values,  and  airspeeds  were  recorded  during  several  maneuvers  performed  with  the 
Tensor  4D  aircraft  with  the  Kestrel  autopilot  installed.  The  data  was  collected  at  a  20  Hz  sample 
rate.  A  total  of  200  samples  were  collected  over  a  10  second  period  for  each  maneuver. 

Commanded  pitch  and  roll  doublets  were  performed.  The  data  was  collected  at  two  wind  tunnel 
airspeeds  with  the  model  mounted  at  the  ball  joint  and  the  CG  at  the  midpoint  of  the 
recommended  aircraft  CG  range.  Airspeeds  were  recorded  during  all  of  the  tests. 

The  pitch  and  roll  angle  measurements  and  servo  displacement  values  measured  by  the  autopilot 
are  provided  as  sample  data.  Plots  of  the  pitch  angle  and  elevator  servo  position  and  roll  angle 
and  aileron  servo  position  are  provided  in  Figs.  5a  and  5b  respectively  for  the  low  velocity  case. 
The  data  trends  recorded  by  the  autopilot  are  consistent  with  the  expected  results.  The  high 
effectiveness  of  the  dual  sets  of  ailerons  on  the  model  can  be  observed  from  the  small  amplitudes 
of  he  roll  aileron  servo  positions.  The  autopilot  is  trying  to  damp  out  the  roll  rate  by  countering 
the  commanded  input  thus  the  servo  position  data  appears  slightly  noisy.  The  pitch  and  roll 
angle  plots  are  smooth  which  is  consistent  with  the  observed  motion.  The  activity  with  the 
vision  based  autopilot  is  mainly  dedicated  to  a  reliable  fast  real-time  connection  between  the 
aircraft  ground  station  located  in  the  wind  tunnel  lab  and  the  vision  camera-large  screen  set 
located  the  Visualization  lab.  The  ground  station  sends  a  set  of  variables  as  the  three  aircraft 
angle  and  rates  and  flight  velocity  to  the  vision  autopilot  that  will  detect  the  motion  through  a 
camera  looking  at  a  virtual  scene  on  a  large  screen.  The  horizon  on  the  screen  moves  with  a 
motion  defined  from  the  real  aircraft  in  the  wind  tunnel.  The  vision  autopilot  reaction  is  sent 
back  to  the  ground  station  through  the  two  laboratories  communication  systems  and  to  the 
aircraft  in  the  wind  tunnel  through  the  wireless  link  between  the  vehicle  and  the  ground  station. 


Figure  5.  a)  Low  speed  pitch  doublet  and  b)  low  speed  roll  doublet.  The  values 
in  radians  are  the  servo  position  (blue)  and  the  aircraft  attitude  (purple). 
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The  originally  desired  frequency  of  data  transfer  was  25  Hz.  This  would  have  allowed  the 
virtual  environment  to  refresh  quickly  enough  to  produce  very  high  quality  imagery  and  to  match 
better  the  expected  dynamics  of  the  aircraft.  However,  a  frequency  constraint  exists  due  to  the 
lack  of  the  GPS  system  on  the  Procerus  autopilot.  The  GPS  system  is  required  for  the  necessary 
packet  to  be  requested  (Packet  18)  in  DevDemo  allowing  the  25  Hz  sampling  rate  transmission. 
Due  to  the  lack  of  GPS,  the  transmission  sampling  rate  was  limited  to  6  Hz. 

An  algorithm  was  implemented  to  create  these  lines  as  well  as  measure  the  pitch  and  roll  from 
the  deviations  and  lines  to  display  the  deviations  from  level  trim  conditions  and  return  them  to 
the  Wind  Tunnel  Laboratory  to  use  for  analysis  [1].  The  camera  sends  the  image  it  ‘sees’  to  the 
workstation  for  the  vision-based  autopilot  to  analyze  as  shown  in  Figure  6.  When  a  yellow  line 
is  visible,  there  is  a  pitch  deviation  from  zero.  When  the  red  line  is  not  horizontal,  there  is  a  roll 
deviation 


Fig.  6.  Virtual  Environment  Display  with  Telemetry  Measurements. 


SAMPLE  OF  RESULTS 

Several  tests  were  performed  in  wind-off  and  wind-on  conditions.  The  output  files  from  the  final 
tests  are  processed  by  a  MATLAB  code  which  creates  a  matrix  from  the  space  delimited  output 
data  file  as  a  time  vector  calculated  dynamically  based  on  the  6Hz  sample  rate.  The  system 
reads  the  aircraft  attitudes  from  the  on-board  IMU  and  estimates  the  same  quantities  from  the 
visual  data  on  the  screen  in  the  visualization  lab.  The  percent  error  for  the  data  is  calculated 
from  the  matrix  of  data  from  the  file  using  the  equation: 
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_  ^  Nominal  -  Experimental 

Percent_Exror  = - - - 100 

Nominal 


Where  the  “Nominal”  value  is  the  telemetry-  data  collected  from  the  aircraft  on-board  system  and 
the  “Experimental”  value  is  the  data  calculated  by  the  vision-based  system.  A  sample  of  results 
is  illustrated  in  Fig.  7  and  Fig.  8  during  a  wind-on  test  with  a  pilot  induced  roll-pitch  doublet 


Aileron  Servo  Position  Data 


Autopilot  Ron 


VisiorvBased  Roll 


Fig.  7.  Wind  Tunnel  Roll  and  Pitch  Test  -  Aileron  Angles  and 
Roll  Data. 
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Elevator  Servo  Position  Data 


Autopilot  Pitch 


Pitrh 


Time,  (s) 


Fig.  8.  Wind  Tunnel  Roll  and  Pitch  Test  -  Elevator 
Angles  and  Pitch  Data. 


CONCLUSIONS 

The  communications  system  created  during  this  phase  of  Task  1  has  resulted  in  a  software 
package  capable  of  testing  a  camera  in  loop  with  an  UAV  operating  in  a  wind  tunnel.  A  series  of 
tests  results  show  that  the  data  processed  in  this  communications  system  is  reliable  and  capable 
of  real  time  operations.  The  maximum  error  between  the  autopilot  telemetry  data  and  data 
calculated  and  sent  from  the  visualization  system  was  less  than  0.3%  with  an  average  value  of 
less  than  0.1%.  The  software  development  package  used  was  Microsoft  Visual  C++  and  resulted 
in  a  GUI  interface,  which  allows  the  user  to  interact  dynamically  with  the  system. 

A  follow-up  phase  could  increase  the  degrees  of  freedom  of  the  aircraft  to  eventually  reach  a  free 
flight  situation  in  the  wind  tunnel  including  an  active  propeller  and  adding  flow  perturbations. 
Finally  a  close  the  loop  control  would  be  implemented  involving  a  feedback  signal  from  the 
visual  control  camera  and  the  UAV  on-board  control  system. 
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Task  2:  Integrated  Visual-Servo  Control  and  Path  Planning  for  Unmanned  Air  Vehicles 

MILESTONES/STATUS 

First  9  Months: 

•  Develop  and  demonstrate  optimization  methods  that  refine  the  path  planning  methods  in 
the  presence  of  obstacles  in  the  image-space  flight  path.  These  algorithms  will  be 
developed  based  on  known  and  unknown  topologies.  (Completed) 

•  Complete  development  and  demonstrate  vision  based  state  estimation  strategies  including 
structure  from  motion  (SFM)  and  simultaneous  localization  and  mapping  (SLAM) 
methods.  (Completed) 

•  Develop  and  demonstrate  visual-servo  controllers  that  use  multi-vehicle  (air  to  air,  air  to 
ground)  images  for  collaborative  mission  objectives.  (Completed) 

Months  10-21 

•  Develop  and  demonstrate  visual-servo  controllers  for  MAV  that  are  robust  to  camera 
uncertainty.  (Completed) 

•  Develop  and  demonstrate  quaternion-based  visual-servo  controllers  to  eliminate  rotation 
singularities.  (Completed) 

•  Develop  and  demonstrate  visual-servo  control  techniques  that  account  for  the  partial 
controllability  of  an  MAV  that  may  be  maneuvering  at  nonlinear  operating  points  away 
from  trim.  (Completed) 

•  Develop  and  demonstrate  optimal  path  planners  for  obstacle  avoidance  that  ensures 
targets  remain  in  the  field  of  view.  (Completed) 

Under  the  leadership  of  Dr.  Dixon  University  of  Florida  and  Dr.  Gans  (UF-REEF)  there  was 
significant  research  accomplishments  achieved:  Further  development  and  verification  of  the 
“Daisy-Chaining”  method  for  large-scale  vision-based  estimation  and  visual  servoing  of 
unmanned  air  and  ground  vehicles  (UAV  and  UGV,  respectively).  Development  of  a  novel 
visual  servoing  method  to  move  a  camera  in  order  to  keep  multiple  objects  in  the  field  of  view 
without  computationally  expensive  image  processing  or  tracking  algorithms  was  achieved. 
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DAISY-CHAINING  VISUAL  SERVOING  FOR  LONG  RANGE  MICRO  AIR  VEHICLE 
(MAV)  CONTROL 

A  principle  concern  in  visual  servoing,  particularly  when  involving  aerial  vehicles,  is  the  limited 
field  of  view  of  the  camera.  When  an  airborne  camera  (such  as  mounted  on  micro  air  vehicle, 
high  altitude  plane,  etc.)  is  used  as  the  sole  sensor  to  generate  control  input  for  a  MAV,  there  are 
several  concerns.  First,  there  must  be  a  method  to  relate  the  position  of  the  MAV  to  fixed 
reference  objects  (e.g.  buildings)  in  order  to  generate  estimates  of  the  current  MAV  pose  and 
velocity  and  relate  it  to  desired  MAV  pose  and  velocity.  Furthermore,  as  the  camera  and  MAV 
move,  fixed  objects  will  leave  the  field  of  view,  and  new  reference  objects  must  be  integrated 
into  the  estimation  algorithm.  This  is  a  difficult  problem.  As  seen  in  Figure  9,  for  one  camera, 
one  MAV  and  one  reference  object,  eight  pose  estimates  must  be  solved  and  chained.  When  a 
second  reference  object  in  involved  (i.e.  a  new  reference  object  has  entered  the  field  of  view)  this 
expands  to  sixteen  pose  estimates.  Furthermore,  not  all  targets  may  be  in  the  camera  field  of 
view  at  the  same  time,  requiring  a  robust  estimation  method  for  this  data. 


Reference 

Object 


Figure  9-Illustration  of  Daisy  Chaining  for  MAV  control 

The  Daisy  Chaining  method  can  be  used  to  solve  this  problem.  By  exploiting  multi-view 
geometry  and  structure- from-motion  methods  of  computer  vision  with  nonlinear  control  theory, 
the  UF  has  designed  a  vision-based  control  system  that  can  regulate  a  MAV  to  a  desired 
trajectory  using  vision  data  from  an  airborne  camera  that  accepts  new  reference  objects  as  they 
enter  the  field  of  view,  and  drops  previous  ones  as  they  leave. 

UNDERDETERMINED  VISUAL  SERVOING  FOR  MULTIPLE  TARGET  TRACKING 

The  vast  majority  of  visual  servo  controllers  are  designed  to  regulate  the  robot,  vehicle,  etc.  to  a 
specific  fixed  pose  or  along  a  specific  fixed  trajectory.  However,  there  are  many  tasks  do  not 
require  or  do  not  have  a  fixed  goal  pose  or  trajectory  for  the  robot.  Thus  a  novel  visual  servo 
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control  method  was  explored  that  does  not  regulate  the  robot  to  a  goal,  rather  the  control 
regulates  some  task  functions  and  allows  the  robot  to  move  to  whatever  pose  is  necessary  to 
regulate  the  task  functions.  This  freedom  of  motion  leads  to  this  controller  being  an 
underdetermined  system.  A  further  benefit  of  the  underdetermined  approach  is  the  expensive 
feature  extraction,  tracking,  and  matching  algorithms  are  not  needed  to  identify  and  associate 
targets  in  each  image. 


The  challenging  task  is  tracking  one  or  more  targets  (Fig.  10  and  Fig.  1 1).  Example  scenarios 
are  keeping  a  moving  vehicle  or  a  dispersing  crowd  in  the  field  of  view  of  an  airborne  camera. 
This  can  be  done  by  regulating  the  distribution  (e.g.  mean  and  variance)  of  targets  in  an  image. 
If  mean  target  position  is  regulated  to  the  image  center  and  variance  of  target  positions  is 
regulated  to  less  than  Vi  the  image  width,  at  least  75  %  of  all  targets  must  be  in  the  field  of  view. 
Depending  on  the  distribution,  this  can  increase  up  to  95%.  Additional  task  functions  can  be 
added,  such  as  keeping  a  desired  distance  or  orientation  to  the  targets.  One  task  the  UF  has 
explored  is  maximizing  perceptibility,  which  is  a  measure  of  how  well  the  airborne  camera  can 
quickly  change  the  distribution  of  targets  in  the  image.  Simulations  have  shown  that  this  method 
is  adept  at  keeping  multiple  objects  in  view,  and  that  perceptibility  help  tracking. 


Fig.  10-  Keeping  multiple  moving 
objects  in  the  field  of  view.  For 
simplicity,  pink  targets  are  added  to 
targets  for  recognition. 


Fig.  11-  Keeping  a  single  moving  object 
in  the  field  of  view.  Target  recognition 
is  used  to  find  the  vehicle  without  using 
artificial  pink  targets. 


VISUAL-SERVO  CONTROL  FOR  UNMANNED  AIR  VEHICLES 

1)  Extension  of  a  novel  visual  servoing  method  to  move  a  camera  in  order  to  track  multiple 
teams  of  targets,  while  maintaining  network  communications  between  multiple  moving  MAV’s. 

2)  Extension  of  the  Hardware  In  the  Loop  Simulation  (HILS)  to  include  a  physical  MAV  in  a 
wind  tunnel,  connected  to  a  the  visualization  system. 

Ensuring  Network  Connectivity  of  MAV's  Performing  Video  Reconnaissance: 

Motivated  by  mission  scenarios  and  sensor  restrictions,  operations  may  require  the  collaboration 
of  assets  over  an  ad-hoc  network.  We  have  explored  the  problem  of  balancing  trade-offs  between 
asset/sensor  cone  positioning  to  satisfy  target  tracking  requirements  and  network  requirements 
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such  as  maintaining  network  connectivity.  There  is  a  direct  need  for  methods  that  manage 
trade-offs  between  asset/sensor  cone  positioning  to  satisfy  mission  requirements  and  network 
requirements  to  ensure  effective  collaboration  between  the  assets.  Yet,  literature  that  focuses  on 
such  issues  appears  sparse. 

The  problem  considered  involves  multiple  MAV’s  equipped  with  cameras.  Each  MAV  is 
dedicated  to  monitoring/tracking  a  team  of  moving  targets,  but  the  communication  network 
between  the  MAV’s  must  be  maintained  for  mission  success.  To  address  the  trade-offs  between 
asset  positioning  and  network  connectivity,  a  prioritized  task-function  based  guidance  law  was 
developed  for  a  simple  scenario  containing  three  assets.  One  developed  task-function  maintains  a 
communication  network  by  ensuring  the  distance  between  the  MAV's  does  not  exceed  a  critical 
threshold.  Additional  task-functions  enable  assets  to  keep  targets  of  interest  in  the  image  cone  by 
regulating  image  features  derived  from  the  camera  view.  Simulations  demonstrated  that  such 
tradeoffs  could  be  successfully  achieved. 


EXTENSION  OF  THE  HARDWARE  IN  THE  LOOP  SIMULATION  PLATFORM  FOR 
VISION-BASED  CONTROL  OF  MTCRO  AIR  VEHICLES 

A  desired  outcome  for  Task  2  is  to  develop  a  sophisticated  simulation  test  bed  for  the  vision- 
based  control  of  MAV’s.  This  test  bed  will  ultimately  provide  multiple  stages  of  increasing 
hardware  interaction.  The  first  stage  is  a  virtual  reality  data  base  capable  of  displaying  large 
environments  and  controls  processing.  The  second  stage  is  a  system  of  modular  displays,  digital 
cameras  and  computers  for  image  and  control  processing.  The  third  stage,  recently 
implemented,  incorporates  a  wind  tunnel,  MAV  and  force  measurement  system.  The  MAV  is 
mounted  on  a  sting  balance  in  a  low  turbulence  wind  tunnel.  Actuation  of  airfoils  in  the  airflow 
of  a  wind  tunnel  will  change  the  attitude  of  the  MAV.  The  attitude  change  is  relayed  to  the  VE 
software,  which  alters  the  viewpoint  accordingly.  In  turn,  the  vision-based  control  routines  send 
velocity  commands  to  the  MAV  autopilot  (see  details  under  Task  1). 

Task  3:  Aerodynamic  and  Structural  Characterization  of  Flexible  and  Morphing  MAV 

MILESTONES/STATUS 

First  9  Months/Status: 

•  Complete  acquisition  of  and  document  experimental  data  sets  (force/moment,  flow  field, 
deformation/shape)  on  fixed  and  passively  flexible  wings.  (Completed) 

•  Collaborate  with  computational  group  (Task  5)  to  assess  and  evaluate  results  of  fixed  and 
passively  flexible  wing.  (Completed) 

•  Define  (with  input  from  Task  4)  and  build  actively  morphing  wing  for  investigating  time- 
dependent  phenomena.  (Completed) 

•  Initiate  acquisition  of  experimental  data  sets  (force/moment,  flow  field,  deformation)  on 
actively  morphed  wings. 


Months  10-21 
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•  Complete  acquisition  of  and  document  experimental  data  sets  (force/moment,  flow  field, 
deformation)  on  actively  morphed  wings.  (Completed  at  Gainesville) 

•  Collaborate  with  computational  group  (Task  5)  to  assess  and  evaluate  results  of  morphed 
wings.  (Collaborated  with  controls  group  on  stability  and  control) 

•  Define  (with  input  from  Task  4)  and  build  actively  flapping  wing  for  investigating  time- 
dependent  phenomena.  (Completed) 

•  Initiate  acquisition  of  experimental  data  sets  (force/moment,  flow  field,  deformation)  on 
actively  flapping  wings.  (Completed) 

GENERAL 

Under  the  leadership  of  Dr.  Ukeiley  and  Dr.  Albertani  from  the  UF-REEF  and  Dr.  Hubner  from 
the  University  of  Alabama  (UA)  there  has  been  significant  progress  in  the  experimental  effort  to 
characterize  the  aerodynamic  and  structural  characteristics  of  MAV  style  airfoils  and  wings. 
This  progress  has  been  on  several  fronts  by  the  team  at  the  UF-REEF  and  the  University  of 
Alabama  which  include; 

•  The  study  of  the  effects  of  passive  flexibility  on  separation  characteristics  of  2-D  airfoils 

•  The  characterization  of  steady  and  unsteady  deformation  of  flexible  wings  including 
modeling  of  aeroelastic  effects 

•  Steady  and  unsteady  aerodynamic  characterization  of  MAV  configurations 

•  Detailed  flow  measurements  of  the  of  low  aspect  ratio  wings 


CHARACTERIZATION  OF  LOW  REYNOLDS  NUMBER  WIND  TUNNEL 

The  Aerodynamic  Characterization  Facility  was  designed  and  installed  at  the  UF-REEF  under 
efforts  partially  supported  by  the  R1AGPS.  This  has  included  a  large  effort  to  quantify  the 
facility  characteristics  which  is  being  reported  in  Albertani  et  al  (2009).  Photographs  of  the 
facility  are  shown  in  12  and  13  below.  The  facility  is  housed  in  a  large  laboratory  of 
approximately  50  x  60  feet  to  provide  a  large  enough  environment  to  isolate  the  effects  of  the 
surroundings  from  the  flow  through  the  wind  tunnel.  As  can  be  seen  in  the  figures  this  is  an 
open  return  wind  tunnel  with  an  open  jet  test  section.  The  air  for  the  facility  is  moved  by  a  50 
horse  power  60  inch  diameter  fan.  The  fan  is  controlled  by  variable  frequency  drive  electronics 
which  allows  for  the  varying  of  the  wind  tunnel  velocity.  The  flow  path  starts  by  entering  the 
conditioning  section  which  has  a  cross  sectional  area  of  10  foot  by  10  foot  through  a  bell  mouth 
inlet.  The  flow  conditioning  section  consists  of  a  radiused  inlet,  metal  honeycomb,  and  then 
several  screen  packs  to  break  up  any  coherent  large  scale  motions.  This  section  is  then  followed 
by  a  4.5  foot  long  straight  settling  region  which  allows  for  all  any  small  scale  organized  motions 
generated  by  the  flow  conditioning  section  to  relax.  This  settling  chamber  is  then  followed  by  a 
square  contraction  section  with  an  area  ratio  of  8  to  1  whose  contour  was  designed  using  the 
tools  discussed  in  Mathew  (2006).  The  exit  of  the  contraction  is  a  42”  square  which  is  the  start 
of  the  open  jet  test  section  and  is  displayed  in  13  (b).  This  open  jet  test  section  is  10  feet  long 


and  has  a  rigid  enclosure  that  can  be  seen  in  12.  The  enclosure  consists  of  two  parts,  an  outer 
shell  which  consist  of  foam  core  wall  office  structure  and  an  inner  structure  made  of  up  of  80/20 
structural  material.  The  overall  dimensions  of  the  enclosure  are  1 1  ft  high,  12  ft  wide  and  15  ft 
long.  At  the  downstream  end  of  the  test  section  is  a  bell  mouth  fiberglass  piece  leading  into  45.6 
inch  square  opening  to  the  diffuser  section.  The  diffuser  section  transitions  from  a  square  cross- 
section  to  a  60  inch  diameter  circle,  over  a  99  inch  axial  distance,  which  is  connected  to  the  fan 
through  a  flexible  coupling.  The  axial  fan  was  manufactured  by  Howden  Buffalo  and  is  spun  by 
a  50  HP  Reliance  Electric  motor.  The  motor  is  controlled  by  a  Toshiba  VF-AS1  high 
performance  inverter.  This  inverter  (Variable  Frequency  Drive,  VFD)  controls  the  frequency  of 
the  electricity  to  the  motor,  which  in  turn  controls  the  rotation  rate  of  the  fan  and  hence  the  flow 
velocity  in  the  test  section.  At  the  end  of  the  motor  housing  is  inline  flow  silencer. 


Figure  12:  Photographs  of  Aerodynamic  Characterization  Facility. 
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A  b 


Figure  13:  Photographs  of  Test  Section’s  Inlet  (a)  and  Exit  (b). 

Assessment  of  Flow  Quality 

A  large  experimental  effort  was  conducted  to  determine  the  properties  of  the  flow  in  the 
Aerodynamic  Characterization  Facility  that  has  included  both  qualitative  and  quantitative 
assessments. 

Since  the  input  to  set  the  rotation  rate  of  the  fan  is  the  frequency  on  the  Variable  Frequency 
Drive  (VFD),  experiments  to  calibrate  this  value  versus  the  air  velocity  at  the  inlet  to  the  test 
section  were  performed.  Several  experiments  have  been  performed  with  a  Pitot-Static  pressure 
probe  situated  in  the  center  of  the  exit  plane  of  the  end  of  the  contraction.  All  of  these 
experiments  involved  setting  the  VFD  to  the  prescribed  frequency  then  letting  the  pressures  in 
tubes  equilibrate  before  taking  the  pressure  data  and  averaging.  This  resulted  in  sampling  the 
pressure  at  2  Hz  for  60  seconds  to  obtain  the  average  velocity  at  the  center  of  the  entrance  to  the 
test  section.  It  should  be  noted  that  all  of  data  presented  here  was  for  VFD  frequencies  of  up  to 
50  Hz  however,  one  can  easily  extrapolate  to  the  maximum  value  of  60  Hz. 
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Figure  14:  Velocity  vs.  Frequency  Calibration  Curve 

Figure  14  displays  a  calibration  curve  to  map  from  the  VFD  to  flow  velocity  at  the  entrance  to 
the  test  section.  These  experiments  were  performed  several  times  with  different  transducers. 
However  the  most  accurate  one  available  at  the  UF-  REEF  is  a  Heise  model  ST-2H  with  both  0.5 
and  2  inches  of  water  transducers  which  have  an  accuracy  of  0.005%  of  the  full  scale  pressure 
reading.  This  results  in  a  velocity  accuracy  of  0.2  m/s  which  are  plotted  as  the  error  bars  on  this 
figure.  Several  experiments  revealed  that  the  measurements  were  repeatable  and  did  not  exhibit 
any  hysteresis  when  the  velocity  was  being  approached  from  either  ascending  or  descending 
values. 

The  Reynolds  number  per  unit  length  and  one  based  on  the  test  section  diameter  are  displayed  in 
Table  1.  This  demonstrates  that  for  an  airfoil  with  a  6  inch  chord  length,  the  Reynolds  numbers 
would  range  from  approximately  16,000  through  175,000.  As  will  be  shown  below  the  facility 
has  been  demonstrated  to  maintain  a  steady  velocity  in  the  1.5  m/s  free  stream  range  validating 
the  overall  goal  of  a  low  Reynolds  number  facility. 


VFD  Setting  (Hz) 

Velocity  (m/s) 

Re  (based  in  test 
section  diameter) 

Re  per  unit 

length  (m'1) 

5 

1.6 

114,000 

107,406 

10 

3.4 

238,000 

220,952 

15 

5.2 

364,000 

336,760 

20 

6.8 

490,000 

454,592 

25 

8.8 

616,000 

574,112 

30 

10.4 

742,000 

692,406 

35 

12.5 

868,000 

810,838 

19 


40 

14.2 

994,000 

929,567 

45 

15.8 

1,120,000 

1,047,037 

50 

17.6 

1,246,000 

1,162,792 

Table  1:  Flow  Velocity  and  Reynolds  Number  versus  VFD  Setting 


Flow  uniformity  studies  were  conducted  with  a  total  pressure  rake  which  had  32  ports  each 
separated  by  1  inch.  The  rake  itself  along  with  the  way  it  was  mounted  is  displayed  in  Figure  15. 
For  the  flow  uniformity  studies  the  rake  was  traversed  (manually)  through  21  vertical  locations 
and  2  spanwise  locations  to  develop  a  uniform  grid  spaced  by  2  inched  in  both  directions.  The 
pressure  readings  were  acquired  with  the  Pressure  Systems  Inc.  scanning  pressure  system  with 
16  ports  with  a  1  psi  full  scale  range.  This  system  quotes  an  accuracy  of  0.05%  of  full  scale 
which  resulted  in  an  accuracy  1.438  m/s.  For  all  of  the  measurements  the  data  was  averaged 
over  120  samples  which  were  acquired  at  a  rate  of  2  Hz.  All  of  the  velocities  presented  in  this 
section  are  normalized  by  the  velocity  at  the  center  location.  It  should  be  noted  that  the  rake 
measurements  will  only  be  reported  for  one  velocity  setting  (approximately  15  m/s)  although 
flow  uniformity  measurements  were  additionally  acquired  at  a  free  stream  velocity  of  2  m/s  and 
are  documented  in  Albertani  et  al  (2009).  Since  the  accuracy  and  range  of  the  Pressure  Systems 
scanning  pressure  system  does  not  allow  for  reliable  interpretation  of  the  data  at  lower  velocities 
a  Scanni-Valve  pressure  switch  was  used  with  the  Heisse  transducer  discussed  above. 


Figure  15:  Total  Pressure  Rake  for  Flow  Uniformity  Studies. 

Contour  maps  of  the  velocity  field  calculated  from  the  total  pressure  measurements  at  4  axial 
locations  are  displayed  in  Figure  16  for  a  40  Hz  setting  on  the  VFD  which  corresponds  to 
centerline  velocity  of  approximately  1 5  m/s.  From  these  plots  the  flow  possesses  a  uniform  core 
that  covers  at  least  60%  of  the  cross  sectional  area  even  at  62.5  %  of  the  axial  extent  of  the  test 
section.  Qualitatively  there  are  some  small  asymmetries  but  in  general  the  flow  appears  to 
evolve  in  a  symmetric  fashion. 

In  order  to  gain  a  better  quantitative  assessment  of  the  flow  quality,  plots  of  the  normalized 
velocities  for  both  centerline  horizontal  and  vertical  slices  at  the  4  axial  positions  are  presented 
in  Figure  17.  From  these  plots  it  is  clear  that  there  are  no  significant  trends  from  top  to  bottom 
or  left  to  right  in  the  wind  tunnel.  Also  apparent  in  these  plots  is  how  small  the  variations  are  on 
the  order  of  less  than  0.5%. 


20 


Y/H 


c)  40  %  of  Test  Section  Length 


Figure  16:  Flow  Uniformity  at  40  Hz  VFD  Setting  (~15  m/s). 


(a)  Vertical 


(b)  Horizontal 


Figure  17:  Centerline  Vertical  and  Horizontal  Velocity  Profiles  at  40  Hz  on  the  VFD  (~15 
m/s). 

Although  the  pressure  measurements  were  quite  encouraging  in  terms  of  flow  uniformity,  a 
device  with  a  higher  frequency  response  was  needed  in  order  to  evaluate  the  unsteadiness.  For 
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the  rest  of  this  section  the  measurements  that  will  be  discussed  obtained  using  a  single 
component  constant  temperature  hot-wire  anemometer  system.  The  sensing  element  of  the  hot 
wire  had  a  diameter  of  5  microns  and  a  length  of  nominally  1  mm.  The  sensing  element  was 
connected  to  a  Dantec  55M10  bridge.  The  output  of  the  bridge  was  then  digitized  with  the 
National  Instruments  4472  card  with  a  sigma-delta  analog-to-digital  converter  with  24  bit 
resolution.  For  these  experiments,  256  blocks  of  2048  points  were  acquired  at  a  rate  of 2048  Hz. 

Table  2  lists  several  quantities  that  were  measured  with  the  hot-wire  anemometer  as  a  function  of 
the  variable  frequency  drive  system  as  displayed  in  the  first  column.  The  second  column  is  the 
mean  velocity  and  the  third  is  standard  deviation  of  the  fluctuating  velocity.  To  calculate 
turbulence  intensity  in  a  wind  tunnel  it  is  common  practice  to  high  pass  filter  the  fluctuations  so 
that  those  associated  with  scales  larger  than  the  facility  are  not  taken  into  account.  After  filtering 
out  these  scales  the  values  are  listed  in  the  last  column  of  Table  2  and  are  all  less  than  0.25%  for 
free  stream  velocities  greater  than  1  m/s. 


VFD 

Setting 

U,  Mean  Velocity 
(m/s) 

<u  >(m/s) 

No  Filtering 

Turbulence 

Intensity 

filtcred>/U 

1 

0.1362 

0.005679 

4.0  % 

3 

0.2895 

0.004152 

1.4% 

5 

0.9492 

0.001213 

0.12% 

7 

2.3154 

0.002255 

0.09  % 

10 

3.3505 

0.004158 

0.06  % 

15 

5.1066 

0.008787 

0.13% 

20 

6.8934 

0.012275 

0.11  % 

25 

8.7058 

0.017555 

0.16% 

30 

10.4996 

0.027106 

0.21% 

35 

12.2955 

0.027439 

0.17% 

40 

14.0959 

0.024822 

0.14  % 

45 

15.8772 

0.024201 

0.07  % 

50 

17.6325 

0.024383 

0.06  % 

Table  2:  Mean  Velocity  and  Turbulence  Intensities 
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Figure  18  displays  the  velocity  power  spectral  density  plotted  at  several  settings  of  the  VFD. 
The  general  trend  of  the  data  from  5  Hz  on  the  VFD  and  up  is  expected,  i.e.,  increasing  energy 
levels  in  the  fluctuating  velocity  as  the  speed  of  the  tunnel  is  increased.  From  this  plot  one  can 
also  view  the  extension  of  the  smaller  scale  (higher  frequency)  phenomena  as  one  would  expect 
with  increasing  Reynolds  number. 


Figure  18:  Power  Spectral  Density  vs.  Frequency  Plots. 

Although  the  overall  behavior  for  the  higher  frequency  cases  shows  the  correct  trends  there  are 
several  peaks  that  might  be  considered  troublesome.  The  first  are  those  that  occur  in  the  range 
of  flow  frequencies  between  3  and  8  Hz  which  begin  to  become  apparent  between  drive 
frequencies  of  1 5  and  20  Hz  and  persist  to  the  maximum  velocities  tested.  To  investigate  these 
peaks  the  power  spectral  densities  are  plotted  versus  Strouhal  number  (fD/U)  in  Figure  19  with 
the  length  scale,  D,  being  the  42  in  exit  length  of  the  contraction.  Clearly,  the  peaks  in  this 
frequency  range  appear  to  have  collapsed  around  a  Strouhal  number  of  0.45.  This  number  is  in 
the  range  of  what  is  generally  expected  for  the  Jet  Column  Mode. 


Figure  19:  Power  Spectra  Plotted  vs.  Strouhal  Number. 

These  measurements  discussed  above  have  demonstrated  the  quality  of  the  flow  in  the 
Aerodynamic  Characterization  Facility  and  it  will  surely  be  a  valuable  facility  for  applied  and 
fundamental  studies  of  the  MAY  aerodynamics.  Additional  measurements  reported  in  Albcrtani 
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et  al  (2008)  along  with  ongoing  further  measurements  demonstrate  that  in  addition  to  steady  flow 
the  facility  will  be  able  to  simulate  gusting  phenomena  that  one  might  find  in  the  urban 
environment. 

AERODYNAMIC  AND  STRUCTURAL  RESEARCH  ACTIVITIES 

Figure  20  outlines  the  aerodynamic  and  structural  research  areas  considered  for  MAV  with 
flexible  and  morphing  wing  structures 

Figure  20.  Aerody  namic  and  Structural  Characterization  of  Flexible  and  Morphing  MAV’s 
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EXPERIMENTS  ON  MICRO-FLAPS 

The  experiments  were  conducted  at  the  low-speed  low-turbulence  wind  tunnel  at  the  University 
of  Florida.  Two  flaps  configuration  with  the  same  height  but  different  spanwise  lengths  were 
tested  on  a  wing  with  a  MAC  of  0.1073  m,  a  wingspan  of  0.377  m,  and  an  aspect  ratio  of  3.78. 
The  airfoil  has  a  mildly  reflexed,  very-thin  section  with  a  thickness  of  0.90  mm  and  a  camber  of 
6  %  at  25  %  MAC.  The  wind  tunnel  tests  consisted  of  angle-of-attack  (AOA)  sweeps  at  73,300 
and  95,200  Reynolds  numbers  based  on  the  MAC  including  the  clean  wing. 

The  experiments  conducted  provided  evidence  that  by  using  a  Gurney  flap  on  a  MAV  wing  there 
exists  a  beneficial  effect  on  lift,  a  reduction  in  drag  in  certain  conditions,  and  an  increase  in 
maximum  lift-to-drag  ratio  compared  to  the  clean  airfoil.  Figure  21  shows  the  increase  of  lift 
coefficient  of  the  wing  with  small  (50%  b)  and  large  (75%  b)  Gurney  flaps  in  respect  to  the  clean 
wing,  at  a  free  stream  Reynolds  number  of  95,200.  The  same  figure  suggests  a  slight  increase  of 
the  lift  curve  slope  for  the  wing  with  flap  and  a  decrease  of  the  stall  angle  and  the  angle  for  zero 
lift  for  the  wings  with  flaps. 


Fig.  21.  Lift  for  clean  wing  and  with 
2.8%  MAC  Gurney  flaps. 


PASSIVE  AEROELASTIC  WINGS  IN  UNSTEADY  CONDITIONS 

Micro  air  vehicles  rigid  and  passive  compliant  wings  were  tested  in  unsteady  conditions  at  the 
REEF  low-speed  low-turbulence  wind  tunnel  equipped  with  a  two  degrees-of-freedom  dynamic 
test  rig.  The  motion,  estimated  by  an  inverse  kinematics  model  that  solves  for  the  desired  testing 
conditions,  is  provided  by  two  ironless  magnetic  linear  motors.  This  facility  permits,  for  the  first 
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time  in  applications  to  micro  air  vehicles,  the  measurement  of  the  two  individual  components  of 
the  rotary  damping  moment.  The  coefficients  due  to  rate  of  change  of  angle  of  attack  were 
measured  by  oscillating  the  model  in  pure  plunging  motion.  The  coefficients  due  to  pitch  rate 
were  measured  during  pitching  motion,  as  illustrated  in  Fig.  22.  To  elucidate  the  correlation 
between  the  wings  structural  deformations  and  the  aerodynamic  damping  characteristics,  rigid 
and  flexible  identical  wings  at  different  levels  of  wing  membrane  tension  were  tested. 


Fig.  22.  Combined  plot  with  the  model  AOA,  pitch  angle  and  lift 
coefficient  versus  time.  Left:  pure  plunging  test.  Right:  pure 
pitching  case.  The  position  of  the  model  is  also  displayed  (light  blue 
line).  The  rectangle  shows  the  area  used  for  data  post-processing. 


The  latex  elastic  membrane  wing  skin  displacements  are  measured  using  the  visual  image 
correlation  therefore  the  pre-tension  strain  state  was  characterized  prior  the  aerodynamic  tests,  as 
illustrated  in  Fig.  23. 
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Fig.  23.  Contour  plots  showing  the  experimental  values  of  the  wings  clastic 
membrane  strain  state  in  wind-off  conditions  for  the  high  pre-tension  case. 

The  strains  are:  a)  in  the  x  direction  and  b)  in  the  y  direction;  c)  and  d)  in  the  low 

The  data  stored  during  the  tests  include  the  time  histories  signals  from  the  six  channels  of  the 
wind  tunnel  sting  balance  and  the  position  of  the  two  linear  motors.  Figure  24  shows  the  lift 
coefficient  time  histories  for  the  case  of  pure  pitching.  The  kinematics  test  conditions  are 
characterized  by  the  pitch  angle  varying  linearly  from  23.01  to  13.0°,  a  constant  angle-of-attack 
(AOA)  of  18.62°  (nominally  18°)  and  a  wind  tunnel  free  stream  velocity  of  10  m/s. 


Time  [ms] 


Time  [ms] 


Fig.  24.  Lift  coefficient  time  histories.  The  plots  represent  the  pure 
pitching  case  (constant  AOA)  for  the  rigid  wing  (left)  and  flexible  wing 
(right).  The  flexible  wing  membrane  is  the  high  pre-tension  condition. 
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The  dynamic  lift  coefficients  compare  properly  with  the  static  counterparts  with  few  anomalies 
as  a  lack  of  lift  when  in  pure-pitching  conditions  at  a  constant  AOA  of  18' . 


UNSTEADY  AERODYNAMICS  AND  STRUCTURAL  DYNAMICS  FOR  FLEXIBLE 
WINGS 

•  In  order  to  facilitate  the  testing  of  the  effects  of  unsteady  features  of  the  MAV  type  airfoils  a 
two-degrees-of  freedom  rig.  Fig.  25  and  26,  was  designed.  This  rig  was  designed  to  facilitate 
the  study  of  the  effects  of  dynamic  rotation,  plunging,  pure  pitching  and  flapping.  These  rigs 
have  been  and  will  continue  to  be  used  in  the  Low  Reynolds  number  wind  tunnel  and  other 
facilities  at  the  REEF  to  further  the  study  of  the  unsteady  flow  effects  on  the  MAV  designs. 


Fig.  25.  The  two-degrees-of  freedom 
dynamic  test  rig  assembled  at  the  side  of 
the  new  wind  tunnel  test  enclosure. 


Fig.  26.  The  one-degree-of  freedom 
dynamic  test  rig  in  operation  in  the  small 
wind  tunnel. 


Wind  tunnel  experiments  on  a  powered  micro  air  vehicle  (MAV)  with  various  flexible  wing 
designs  and  variable  elevator  deflections  proved  effective  for  obtaining  the  aerodynamic 
coefficients  in  forms  of  linear  regression  models  and  investigating  the  mutual  dependencies 
between  the  variables  [2]. 

Models  for  the  aerodynamic  and  propulsion  coefficients,  plotted  in  Figure  27,  were  obtained 
from  multi-factorial  wind  tunnel  experiments  the  nonlinear  dependences  of  the  propeller  speed 
on  the  angle-of-attack,  free  stream  velocity  and  elevator  angle. 
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Fig.  27.  Aerodynamic  coefficient  response  surfaces  for  6, 7,  and  8  volts  of  input  power  to  the  DC 
micro  motor. 
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Figure  28.  Time-averaged  membrane  displacements  (mm)  of  the  wing  with  a  high  pre-tension. 


Fig.  29.  Average  spectrum  (velocity/force)  from  LDV  testing  of  membrane  wing  with  high  pre¬ 
tension. 


Experiments  were  also  carried  out  to  evaluate  the  fundamental  structural  dynamics  parameters 
(mode  shapes,  natural  frequencies,  and  modal  damping)  of  membrane  MAV  wings  [3].  Bench 
vibrations  tests  (dry  tests)  and  naturally  induced  vibrations  measurements  in  a  wind  tunnel  were 
performed  on  15  cm  wings  with  different  levels  of  elastic  membrane  pre-tension.  The  average 
wing  membrane  out-of-plane  deformations  in  wind  tunnel  tests  are  presented  in  Figure  28  for 
three  angles  of  attack  and  at  a  free  stream  velocity  of  9.5  m/s.  The  average  spectrum  in  terms  of 
velocity/force  in  natural  induced  vibrations  is  depicted  in  Figure  29  for  the  case  with  high 
membrane  pre-tension.  Results  show  an  acceptable  experimental  accuracy  for  the  modal 
parameters. 
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An  experimental  method  based  on  dynamic  visual  image  correlation  and  homogeneous 
transformation  matrix  post-processing  was  developed  for  decoupling  the  kinematics  of  wing 
motion  from  the  deformation  of  a  flapping  flexible  wing  [4].  Figure  30  displays  the  combined 
wing  twist  and  elastic  deformation  occurring  at  the  10  Hz  flapping  frequency  of  a  latex 
membrane  flexible  wing. 


Flapping  Frequency:  10Hz  Flapping  Frequency:  10Hz  w,imn,l 


Fig.  30,  Projected  rigid-body-motion  plotted  against  measured  displacement  for  a  membrane 
latex  flapping  wing  at  a)  Start  of  up-stroke  and  b)  Start  of  down-stroke. 

The  primary  objective  of  the  UA  effort  was  to  study  how  a  flexible  (latex)  membrane,  similar  to 
those  used  on  UF-designed  MAVs,  affects  the  wake  structure  of  flat  and  cambered  plates  at  low 
Re.  Past  and  recent  research  has  shown  that  passive  geometric  changes  and  flexibility  in  the 
wing  structure  can  improve  MAV  performance  and  alleviate  sensitivity  to  flow  disturbances.  At 
higher  angles-of-attack  (AOA),  the  extension  and/or  passive  vibration  of  the  membrane  could 
assist  in  the  reattachment  of  the  leading-edge  shear  layer  and  decrease  wake  width  and  drag. 
Figures  31  and  32  shows  the  two  basic  wing  designs  studied:  perimeter  reinforced  (PR)  and 
batten  reinforced  (BR).  The  latter  does  not  have  a  rigid  trailing-edge  and,  thus,  does  not  prevent 
vibration  of  the  membrane  at  the  trailing-edge.  High  aspect  ratio  plates  (AR  =  8),  designed  to 
minimize  the  tip  effects,  with  multiple  membrane  cells  and  rounded  leading-edges  were  tested  at 
Re=45,000. 

Surface  and  flow  visualization  confirmed  leading-edge  separation  and  reattachment  on  the  solid 
flat  plate  for  small  AOA  (see  figure  32).  The  separation  bubble  increased  with  AOA,  becoming 
fully  separated  at  angles  greater  than  8°.  For  the  cambered  plates,  the  onset  of  leading-edge 
separation  occurred  at  a  positive  AOA  but  slightly  negative  local  leading-edge  AOA  relative  to 
the  free  stream  flow.  Based  on  the  jump  in  wake  width,  the  flow  appeared  fully  separated  on  the 
cambered  plate  at  6°  AOA.  The  hot-wire  wake  surveys  quantified  velocity  deficit  and  turbulence 


30 


Fig.  31.  Flat  and  cambered  (6%)  membrane  p.g  J2  Separa|io„  bubbIe  on  a  flat  plaK  at  2. 

intensity  profiles  as  well  as  width  and  local  drag  coefficients.  Fig.  33  is  a  comparison  of  the  flat 
plate  wake  width  at  50%  velocity  deficit  (left)  and  local  drag  coefficient  (right)  for  the  flat  plate. 
For  the  PR  wing,  there  was  a  slight  increase  in  wake  width,  local  drag  and  turbulence  intensity 
for  AOA  <  6°.  Slight  improvement  occurred  at  7°,  but  the  global  characteristics  were  relatively 
the  same  at  8°  (near  full  separation).  The  wake  structure  was  vastly  different  for  the  BR  wing  due 
to  trailing-edge  flapping.  While  the  wake  deficit  was  lower  than  that  of  the  PR  and  solid  wings, 
the  wake  width,  local  drag  and  turbulence  intensities  were  significantly  increased,  especially  at 
small  AOA.  By  8°,  the  global  characteristics  were  similar  to  that  of  the  PR  and  solid  plates. 
Scalloping  the  trailing  20%  of  the  BR  membrane  showed  considerable  improvement  in  the  local 
drag  coefficient  (based  on  the  shortened  local  chord).  This  was  due  primarily  to  the  substantial 
decrease  in  the  wake  deficit  magnitude.  The  wake  width,  however  was  greater  than  that  of  the 
solid  and  PR  wing,  indicating  trailing-edge  fluctuation  was  still  present  but  not  as  large  as  the 
unscalloped  BR  membrane.  For  the  cambered  plate,  the  presence  of  the  PR  and  BR  membrane 
decreased  the  wake  width  and  local  drag  coefficient  at  6°  AOA  (near  full  separation),  but  showed 
no  improvement  at  lower  or  higher  AOAs.  Similar  to  the  flat  plate,  scalloping  the  BR  membrane 
on  the  camber  plate  significantly  decreased  the  width  and  drag  over  the  pre-  and  post-fully 
separated  flow  range  (2°-14°). 


Figure  33.  EfTect  of  membrane  geometry  on  wake  width  (left)  and  local  drag  coefficient  (right)  for  AOA  - 

solid,  PR,  BR,  BR-scalloped  flat  plate. 
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MAY  WING  VELOCITY  FIELD 


A  series  of  detailed  flow  measurements  over  a  rigid  low  aspect  ratio  MAV  wing  were  conducted 
using  the  REEF’s  time  resolved  particle  Image  Velocimetry  system  and  low  speed  wind  tunnels. 
The  experimental  effort  has  been  reported  in  Khambatta  (2008)  and  only  a  small  set  of  the  data 
will  be  discussed  below.  In  addition  we  have  used  these  results  to  compare  to  some  steady  CFD 
results  and  was  reported  in  Khambatta  et  al  (2008).  The  wing  chosen  for  these  studies  was  the 
rigid  counterparts  to  those  used  for  the  wing  deformation  studies  mentioned  elsewhere  and  were 
developed  by  P.  Ifju  and  his  students.  Below  we  will  detail  the  comparison  between  the  P1V  and 
CFD  results  after  a  brief  review  of  some  details  from  the  two  techniques.  Readers  interested  in  a 
more  elaborate  discussion  of  the  detailed  PIV  measurements  should  review  Khambatta  2008. 

The  figure  34  below  shows  the  planes  in  which  the  PIV  measurements  were  acquired. 
Measurements  for  both  configurations  were  taken  at  a  free-stream  velocity  of  9.5  m/s 
corresponding  to  a  Reynolds  number  of  75,000  based  on  the  wing  chord  length  of  0.124  m. 
Free-stream  conditions  were  verified  before  the  experiment  by  means  of  a  pitot-static  probe 
connected  to  a  Heise  (ST-2H)  0-0.5”  of  H2O,  pressure  transducer  with  an  accuracy  of  0.015%  FS 
(~0.1  m/s).  The  time  scale  to  ensure  the  capture  of  statistically  independent  data  was  calculated 
to  be  ~0.028  s,  assuming  that  the  largest  conceivable  scale  is  on  the  order  of  the  wing-chord. 
This  implies  that  the  sampling  frequency  should  be  less  than  35.5  Hz,  to  ensure  statistically 
independent  samples. 


Figure  34:  Measurement  planes  for  PIV 
experiments 

The  instantaneous  measurements  of  the  velocity  field  were  acquired  using  a  Dantec  Dynamics 
high  speed  PIV  system.  The  PIV  components  include  a  Series  800-PIV/40G  Nd:YAG  laser 
system,  two  IDTX S-5  high  speed  digital  cameras  with  f-2.8  Nikkor  lenses  and  version  4.71  of 
the  Flow  Manager  software.  The  cameras  have  pixel  dimensions  of  12pm  x  12pm  and  a 
resolution  of  1260x1024.  For  the  2-component  PIV  setup,  the  camera  lens  magnification  was  M 
=  0.11  and  for  the  3-component  PIV  setup,  the  average  lens  magnification  was  M  -  0.09.  A 
Laskin  nozzle  seeding  generator  from  Dantec  Dynamics  (model  10F03)  was  used  to  atomize 
olive  oil  and  produce  seeding  particles  expected  to  be  in  the  size  range  of  1-5  pm.  The  seeding 
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particles  are  introduced  into  the  flow  just  upstream  of  the  inlet-contraction  to  the  wind  tunnel  in 
such  a  manner  that  they  produced  a  region  of  uniformly  seeded  flow  all  around  the  wing  within 
the  test  section.  Given  the  magnification  of  the  lenses,  an  J#  =  2.8,  and  the  wavelength  of  the 
laser,  Xl  -  532  nm,  one  can  calculate  the  diffraction  limited  spot  size,  dsp  for  a  single  lens 
(Meinhart  and  Wereley,  2003)  using  the  following  relationship, 

d,r  *  2.44(A/  +  .  (1) 


The  particle  image  diameter,  dr  on  the  pixel  array  is  then  calculated  using  equation  (2),  where  dp 
is  the  diameter  of  the  seeding  particles  (Westerweel,  2000): 

(2) 

For  these  calculations,  dp  was  assumed  to  be  5.0  microns,  based  on  specifications  provided  by 
Dantec.  Based  on  equations  (1)  &  (2),  dsp  &  dr  were  calculated  to  be  4  microns.  For  cases  where 
M  «  1 ,  dr  — *  dsp  (Westerweel,  2000),  and  therefore,  our  ratio  of  particle  image  diameter  to  pixel 
dimension  was  0.33. 

For  2-component  PIV,  the  width  of  the  image  window  was  0.14  m  and  to  be  able  to  apply  an 
interrogation  area  of  16pix  x  16pix  (3.55  mm  square)  to  an  image  pair,  the  time  delay  between 
successive  images  in  a  pair  was  set  at  35  ps.  For  3-component  PIV,  the  width  of  the  image 
windows  were  0.1  m  and  an  interrogation  area  of  32pix  x  32pix  was  used  (2.54  mm  square).  To 
obtain  statistically  independent  data,  the  frame  rate  was  set  at  35.0  Hz.  A  sample  of  the 
convergence  is  taken  from  the  measurements  at  a  point  in  the  freestream  at  a  =  15°  and  is  shown 
in  Figure  35,  the  data  has  been  normalized  by  the  mean  velocity  at  that  point  (calculated  from 
810  samples).  It  can  be  seen  that  once  the  sample  size  exceeded  600,  the  variations  in  the  mean 
velocity  were  less  than  0.2%  of  the  mean.  After  the  raw  images  had  been  obtained,  a  low  pass 
Gaussian  filter  with  a  width  of  0.1  was  applied  to  the  data  and  the  images  were  then  processed 
with  50%  overlap.  A  series  of  high  speed  measurements  were  acquired  at  400  Hz  in  an  effort  to 
resolve  the  temporal  behavior  of  the  more  energetic  features  of  the  flow.  While  the  system  is 
capable  of  sampling  at  up  to  512  Hz  in  double-frame  mode,  measurements  collected  above  400 
Hz  appeared  to  suffer  from  inadequate  optical  conditions.  These  high  speed  measurements  were 
collected  under  limited  conditions:  at  a  =  25°  &  30°  in  the  x-y  plane  at  2b/z  =  0  &  1  (2- 
component  PIV);  at  a  =  30°  in  the  z-y  plane  at  c/x  =  0.5,  0.87  &  1.1  (3-component  PIV). 
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Figure  35:  Convergence  plot  of  normalized  mean  vs.  sample  size. 

All  measurements  (2-component,  3-component  configurations,  statistically  independent  and  high 
speed),  comprised  810  samples.  These  measurement  locations  described  should  provide  a 
sufficient  description  of  the  flow  topology  surrounding  this  MAV  wing  to  be  compared  with 
numerical  simulations,  and  provide  some  insight  into  the  non-linear  lifting  characteristics  of  low 
Reynolds  number  airfoils,  in  particular,  the  vortices  near  the  wing  tips  and  the  regions  of 
separated  flow  on  top  of  the  wing.  Table  3  and  Table  4  present  the  locations  where  PIV 
measurements  were  taken  around  the  wing. 


|  Angles  of  attack  |  0°  |  10°  |  15°~T  18°  I  20°  |  21°  1  25°  |  30°  |  |  1 


Span-wise  locations,  2z/b 
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0.33 

0.5 

0.67 
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0.92 

1.0 

1.08 

plane) 

Table  3:  Angles  of  attack  and  span-wise  locations  for  2-component  PIV  measurements. 
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0.4 
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1.4 

Table  4:  Angles  of  attack  and  chord-wise  locations  for  stereoscopic  PIV  measurements. 
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NUMERICAL  METHODS 


The  three-dimensional,  incompressible  Navier-Stokes  equations  written  in  curvilinear 
coordinates  are  solved  for  the  steady,  laminar  flow  over  the  MAV  wing.  The  fuselage, 
stabilizers,  and  propeller  are  not  taken  into  account.  The  computational  domain  can  be  seen  in 
Figure  36,  with  the  MAV  wing  (red)  enclosed  within.  Inlet  and  outlet  boundaries  are  marked  by 
the  flow  vectors:  velocity  is  specified  at  the  inlet,  while  a  zero  pressure  condition  is  enforced  at 
the  outlet.  The  configuration  shown  in  Figure  36  is  for  simulations  at  a  model  inclination  of  0°. 
For  non-zero  angles,  the  lower  and  upper  surfaces  will  also  see  a  mass  flux.  The  side  walls  are 
modeled  as  slip  walls,  and  thus  no  boundary  layer  forms.  The  dimensions  of  the  computational 
domain  are  given  in  terms  of  the  root  chord,  and  are  placed  far  enough  away  from  the  MAV 
surface  so  as  not  to  significantly  affect  the  aerodynamics.  As  no  flow  is  expected  to  cross  the 
root  chord  of  the  wing  (no  propeller  is  modeled),  symmetry  is  exploited  by  modeling  only  half  of 
the  computational  domain  (the  plane  of  symmetry  is  also  modeled  as  a  slip  wall).  A  detailed 
view  of  the  resulting  structured  mesh  is  given  in  Figure  36.  210,000  nodes  fill  half  of  the 
computational  domain,  with  1300  nodes  on  the  wing  surface. 

In  order  to  handle  the  arbitrarily  shaped  geometries  of  a  micro  air  vehicle  wing,  the  Navier- 
Stokes  equations  must  be  transformed  into  generalized  curvilinear  coordinates:  ^(x,y,z),  q(x,y,z), 
£(x,y,z).  This  transformation  is  achieved  by  (Tannehill  et  al,  1997): 
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where /j  are  metric  terms,  and  J  is  the  determinant  of  the  transformation  matrix: 


,  3(x,y,z) 

3(^,0 


(4) 


Using  the  above  information,  the  steady  Navier-Stokes  equations  can  then  be  written  in  three- 
dimensional  curvilinear  coordinates  (Shyy,  1994).  The  continuity  equation  and  u-momentum 
equation  are  presented  here  in  strong  conservative  form,  with  the  implication  that  the  v-  and  w- 
momentum  equations  can  be  derived  in  a  similar  manner. 


u5+vn+w(.=o 


(5) 
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where  p  is  the  density,  p  is  the  pressure,  p  is  the  viscosity,  q,,  are  parameters  dictated  by  the 
transformation  (expressions  can  be  found  in  Shyy,  1994),  and  U,  V,  and  W  are  the  contravariant 
velocities,  given  by  the  flux  through  a  control  surface  normal  to  the  corresponding  curvilinear 
directions: 


U  =  /,-u+/j-v+^,-w 

V  =  /n-U  +  /n-V  +  /2J-W  (7) 

W=/3.  U+/32-V  +  /33-W 

In  order  to  numerically  solve  the  above  equations,  a  finite  volume  formulation  is  employed, 
using  both  Cartesian  and  contravariant  velocity  components  (Tannehill  et  al,  1997).  The  latter 
can  evaluate  the  flux  at  the  cell  faces  of  the  structured  grid  and  enforce  the  conservation  of  mass. 
A  second  order  central  difference  operator  is  used  for  computations  involving  pressure/diffusive 
terms;  a  second  order  upwind  scheme  handles  all  convective  terms  (Thakur  et  al,  2002). 


Figure  36:  CFD  computational  domain  (left)  and  detail  of  mesh  near  wing  surface  (right). 


Figure  37  through  Figure  41  allow  for  a  qualitative  comparison  and  discussion  about  the 
numerical  and  experimental  results  from  2-component  PIV  at  a  =  10°,  15°,  21°  and  30°  and  at 
2z/b  =  0  &  1.  Due  to  the  dihedral  in  the  wing  and  because  the  camera  in  the  2-component  PIV 
setup  was  perpendicular  to  the  x-y  plane,  the  flow  adjacent  to  the  wing  surface  could  not  be 
characterized  in  certain  regions,  (as  is  indicated  by  a  white  space  between  wing  profile  and  the 
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data).  These  regions  do  not  exceed  0.04c  in  height  above  the  wing  surface.  The  PIV  results 
were  obtained  with  one  digital  camera  capturing  images  at  35  Hz  at  a  resolution  of  1260  x  1024 
pixels  and  at  a  distance  of  0.625  m  from  the  plane  of  illumination. 


(A)  (B) 


Figure  37:  A).  Percentage  distribution  of  estimated  vectors  for  a  =  10°  &  2z/b  =  0,  and  B). 

2z/b  =  1. 

For  a  =  10°  and  2z/b  =  0,  Figure  37  A,  shows  that  the  number  of  estimated  vectors  is  less  than 
10%  of  the  total  number  of  samples  for  the  majority  of  the  area  over  the  wing.  The  data  near  the 
leading  edge  has  a  region  where  the  percentage  of  estimated  vectors  gets  significantly  larger,  but 
this  may  be  attributed  to  the  fact  that  there  is  poor  correlation  between  the  image  pairs  due  to 
those  locations  being  close  to  the  boundary  of  the  image  frame.  Figure  37  B  shows  a  similar 
distribution  of  the  percentage  of  estimated  vectors  for  a  =  10°  and  2z/b  =  1 .  It  is  observed  that 
there  are  a  large  percentage  of  estimated  vectors  in  a  thin  region  that  is  adjacent  to  the  perimeter 
of  the  wing-tip.  This  is  attributed  to  the  edge  of  the  wing-tip  that  reflected  a  portion  of  the  light 
from  the  laser  sheet  illuminating  the  x-y  plane  immediately  adjacent  to  it,  causing  a  bloom  effect 
on  the  sensor  array.  The  estimated  vectors  for  the  other  cases  have  a  similar  distribution  over  the 
vector  map  as  the  cases  shown  in  Figure  37. 

Looking  at  the  time-averaged  data  in  Figure  38  through  Figure  41,  there  is  indication  of  an 
acceleration  of  the  flow  over  the  wing  at  the  point  of  greatest  camber,  as  expected,  after  which 
there  is  some  indication  of  separation  as  the  flow  over  the  wing  comes  in  contact  with  an  adverse 
pressure  gradient.  The  numerical  and  experimental  data  appear  to  correlate  well  on  this  aspect. 
However,  the  numerical  data  indicates  a  stronger  region  of  re-circulation  sitting  within  the 
reflexed  region  of  the  airfoil  (~x/c  =  0.75),  which  is  not  visible  in  the  experimental  data  (for  a  = 
10°  &  15°  at  z/b  =  0).  The  flow  over  the  reflexed  region  is  clearly  slowed  down,  and  this  might 
be  indicative  of  a  smaller  separation  bubble  slightly  more  upstream  than  where  the  numerical 
model  predicts  it  to  be.  Furthermore,  for  a  =  21°,  the  numerical  model  predicts  large  scale 
separation  over  the  wing  and  the  experimental  results  indicate  a  much  smaller  region  of  re¬ 
circulation  sitting  just  above  the  reflexed  portion  of  the  wing.  At  a  =  30°,  both,  the  numerical 
model  and  the  experimental  data  indicate  a  large  scale  separation  over  the  wing,  but  with  varying 
results  for  the  centers  of  the  vortex  cores  at  2z/b  =  0. 
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For  a  =  10°,  15°  &  21°  and  2z/b  =  1,  the  numerical  model  consistently  predicts  a  slower  velocity 
in  the  vortex  tube  as  compared  to  the  experimental  data.  The  flow  in  this  region  is  highly  three 
dimensional,  and  it  is  possible  that  the  experimental  data  was  collected  in  a  region  where  the 
strongest  component  of  the  velocity  was  in  a  stream-wise  direction  rather  than  normal  to  it. 
However,  as  can  be  observed  in  Figure  38  through  Figure  41,  the  stream-wise  component  of  the 
numerical  data  is  also  consistently  slower  in  the  vortex  core  as  compared  to  the  measured  data. 
In  the  cases  where  2z/b  =  1  in  Figure  38  through  Figure  41,  the  stream-wise  velocity  in  the  core 
of  the  separation  region  appears  to  be  higher  than  in  the  free-stream.  This  clearly  indicates  a 
discrepancy  between  the  numerical  simulation  and  the  experimental  results.  However,  the 
experimental  results  were  collected  -2-3  mm  off  the  wing  tip,  (to  avoid  excessive  blooming),  in 
a  region  adjacent  to  the  core  of  the  separated  flow  region,  where  the  stream-w  ise  component  of 
the  velocity  is  greater  than  in  the  core.  In  a  case  where  the  vortex  is  -12  mm  wide  (Figure  42) 
and  the  core  is  only  3-4  mm  wide,  taking  into  account  he  laser  sheet  thickness  (-2  mm),  it  is 
possible  that  the  numerical  and  experimental  data  were  collected  from  two  different  x-y  planes. 


instantaneous  time-averaged  CFD 
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Figure  38:  Normalized  streamwise  flow  at  a  =  10°:  instantaneous  PIV  (left),  time-averaged 
PIV  (center),  and  CFD  (right).  Span  station  2z/b  =  0  (top)  and  2z/b  =  1  (bottom). 
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Figure  39:  Normalized  streamwise  flow  at  o  =  15°:  instantaneous  PIV  (left),  time-averaged 
PIV  (center),  and  CFD  (right).  Span  station  2z/b  =  0  (top)  and  2z/b  =  1  (bottom). 


Figure  40:  Normalized  streamwise  flow  at  a  =  21°:  instantaneous  PIV  (left),  time-averaged 
PIV  (center),  and  CFD  (right).  Span  station  2z/b  =  0  (top)  and  2z/b  =  1  (bottom). 
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Figure  41:  Normalized  stream-wise  flow'  at  a  =  30°:  instantaneous  PIV  (left),  time- 
averaged  PIV  (center),  and  CFD  (right).  Span  station  2z/b  =  0  (top)  and  2z/b  =  1  (bottom). 
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Figure  42:  Normalized  stream-wise  flow  (top)  and  vorticity  (bottom)  at  x/c  =  0.51,  a  =  15°, 
with  normal  components  indicated  by  the  vector  plot:  time-averaged  PIV  (left),  and  CFD 

(right). 
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Figure  43:  Normalized  stream-wise  flow  (top)  and  vorticity  (bottom)  at  x/c  =  1.02,  a  =  15°,  with  normal 
components  indicated  by  the  vector  plot:  time-averaged  PIV  (left),  and  CFD  (right). 
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Figure  44:  Normalized  stream-wise  flow  (top)  and  vorticity  (bottom)  at  x/c  =  0.51,  a  =  21°, 
with  normal  components  indicated  by  the  vector  plot:  time-averaged  PIV  (left),  and  CFD 

(right). 
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Figure  45:  Normalized  stream-wise  flow  (top)  and  vorticity  (bottom)  at  x/c  =  1.02,  a  =  21°, 
with  normal  components  indicated  by  the  vector  plot:  time-averaged  PIV  (left),  and  CFD 

(right). 


Figures  42  through  Figure  45  show  results  from  the  3-component  PIV  tests  and  the 
corresponding  numerical  data.  During  the  collection  of  the  experimental  data,  the  wing  had  been 
installed  with  a  slight  roll  angle  with  respect  to  the  cameras.  In  an  effort  to  correctly  depict  and 
compare  the  flow  around  the  wing  with  a  roll  angle  of  0°,  the  co-ordinate  system  and  the  flow 
field  have  been  rotated,  thereby  giving  the  impression  of  a  skewed  flow  field. 

The  wing  tip  vortices  in  these  cases  can  be  clearly  seen  and  when  observed  over  a  period  of  time, 
appear  to  grow  and  shrink  in  size,  alternating  between  wing-tips.  In  the  a  =  21°  case,  the  region 
of  separation  on  top  and  behind  the  wing  is  clearly  smaller  in  magnitude  in  the  experimental  data 
than  in  the  numerical  data,  and  this  correlates  with  the  comparison  done  earlier  with  the  2- 
component  PIV  data,  where  the  numerical  model  predicted  massive  separation  and  the 
experimental  data  showed  otherwise.  A  quantitative  comparison  has  been  drawn  between  the 
numerical  and  experimental  results  in  Table  5,  based  on  the  locations  of  the  vortex  core  centers 
and  the  vorticity  at  that  location. 


The  downward  shift  in  the  locations  of  the  vortex  cores  with  an  increase  in  z/c  is  due  to  the 
downwash  caused  over  the  wing  due  to  its  angle  of  attack  and  geometry  .  Not  only  is  there  a 
downwash  effect,  there  also  appears  to  be  a  side-ways  shift  of  the  vortex  core  towards  the  center 
of  the  wing,  this  trend  is  observed  in  both,  the  numerical  and  the  experimental  data.  The 
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vorticity  at  the  wing-tips  is  also  expected  to  increase  with  an  increase  in  angle  of  attack, 
however,  the  experimental  data  indicates  otherwise,  showing  a  decay  in  the  vorticity  at  the 
vortex  center  at  a  =  21°  when  compared  at  a  =  15°. 


Angle  of 
attack,  a 

Chord-wise 
location  (z/c) 

Location  of  vortex  core 

Vorticity  (r1) 

Z 

(mm  from  origin) 

y 

(mm  from  origin) 

15° 

0.5 

Experimental 

78.8 

2.1 

4049 

Numencal 

75.5 

4.0 

3112 

15° 

1.0 

Experimental 

75.9 

-8.8 

3593 

Numerical 

71.4 

-8.1 

1313 

21° 

0.5 

Experimental 

79.0 

-0.4 

3579 

Numerical 

75.5 

-1.4 

3681 

21° 

1.0 

Experimental 

74.4 

-13.8 

3508 

Numerical 

71.4 

-19.1 

1332 

Table  5:  Comparison  of  vortex  core  locations  and  vorticity,  between  experimental  and 

numerical  results. 
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Task  4:  Control  of  Biologically-Inspired  Morphing  for  Variable-Geometry’  MAV 
MILESTONES/STATUS: 

First  9  Months/Status: 

•  Complete  analysis  of  aerodynamic  shape-varying  configurations  (collaborate  with  REEF 
wind  tunnel  lab  and  REEF  CFD  researchers).  (Completed) 

•  Analyze  structural  dynamics  of  vehicle  with  morphing  implementations  (Collaborate 
with  AVCAAF  researchers).  (Completed) 

•  Develop  framework  to  represent  shape-dependent  models  (collaborate  with  GNV 
researchers),  (completed) 

•  Conduct  stability  analysis  for  configurations  within  the  framework  (collaborate  with 
GNV  and  AFRL  researchers).  (Completed) 

•  Implementation  of  configurations  to  alter  CG  and  tail  (collaborate  with  REEF  rapid 
prototyping  lab).  (Completed) 


GENERAL 

Under  the  leadership  of  Dr.  Rick  Lind  from  the  University  of  Florida  and  Dr.  Chakravarthy 
from  the  UF-REEF  significant  progress  has  been  made  in  researching  the  effects  of  time- 
varying  inertias  on  flight  dynamics  of  an  asymmetric  variable-sweep  morphing  micro  air 
vehicles.  Mission  capability  can  be  enabled  by  morphing  a  micro  air  vehicles  to  optimize 
its  aerodynamics  and  associated  flight  dynamics  for  each  maneuver.  Such  optimization 
techniques  often  consider  the  steady-state  behavior  of  the  configuration;  however,  the 
transient  behavior  must  also  be  analyzed.  In  particular,  the  time-varying  inertias  have  an 
effect  on  the  flight  dynamics  that  can  adversely  affect  mission  performance  if  not  properly 
compensated.  These  inertia  terms  cause  coupling  between  the  longitudinal  and  lateral- 
directional  dynamics  even  for  maneuvers  around  trim.  A  simulation  of  a  variable-sweep 
aircraft  undergoing  a  symmetric  morphing  for  an  altitude  change  shows  a  noticeable  lateral 
translation  in  the  flight  path  because  of  the  induced  asymmetry.  None  linear  dynamic 
equations  including  varying  mass  distributions  where  derived  and  applied  to  a  variable 
wing  sweep  vehicle  described  below: 

TIME-VARYING  DYNAMICS  OF  MORPHING  AIRCRAFT 

A  study  into  the  time-varying  dynamics  of  the  morphing  wing,  double  section,  variable-sweep 
MAV  has  been  in  progress.  In  this  study,  the  time-varying  counterparts  of  the  conventional 
short  period,  phugoid,  roll,  spiral  and  dutch  roll  modes  were  researched.  More  specifically,  if  the 
MAV  is  hit  by  a  disturbance  at  the  same  time  while  it  is  executing  a  morphing  trajectory,  this 
study  characterizes  the  modes  of  the  MAV  at  such  times. 

It  is  well  known  that  the  study  of  time-varying  systems  (even  linear  ones)  is  far  more  complex 
than  those  of  LTI  systems.  For  a  LTV  system,  since  the  eigenvalues  of  the  system  matrix  in 
general  have  no  significance  with  regard  to  providing  information  of  the  system’s  stability  or  its 
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modal  behavior,  alternative  techniques  were  explored.  Several  techniques  have  been  proposed  in 
the  literature,  but  all  of  them  suffer  from  some  drawbacks  or  the  other.  Wu  (1984)  proposed  a 
system  of  using  (extended)  x-eigenvalues  and  x-eigenvectors  (using  new  definitions  for  these 
concepts);  this  method  is  good  in  that  it  can  always  diagonalize  the  system  matrix  (resulting  in 
decoupled  modes);  however  it  has  the  drawback  that  the  x-eigenvalues  and  x-eigenvectors  are 
not  unique.  Indeed,  each  x-eigenvalue,  x-eigenvector  pair  makes  some  sense  only  as  a  pair,  but 
not  individually.  Also,  if  the  time-variation  is  for  a  finite  duration,  these  x-eigenpairs  do  not 
necessarily  merge  with  the  LTI  eigenvalues  and  eigenvectors  that  occur  at  the  start/end  of  the 
time-variation.  Kamen  (1988)  proposed  a  factorization  approach  to  determine  the  time-varying 
poles  (and  thus  the  modes)  of  an  LTI  system.  This  method  involves  the  solution  of  a  nonlinear 
time-varying  Riccatti  equation.  This  method,  when  it  works,  is  most  elegant;  in  particular,  its 
time-varying  poles  can  be  related  with  the  LTI  poles  that  occur  at  the  start  time  of  a  finite 
duration  morphing,  but  at  other  times,  the  Riccatti  equation  can  have  unstable  solutions,  and  this 
causes  the  time-varying  poles  to  diverge  to  negative  oo.  Zhu  (1991)  proposed  notions  of  parallel 
D  spectra  and  scries  D  spectra,  and  used  these  notions  to  diagonalize  the  system  matrix  (in 
addition  to  transforming  to  other  canonical  forms)  but  they  too  have  the  problem  that  the  time- 
varying  poles  can,  under  certain  conditions  diverge  to  negative  oo.  O’Brien  and  Iglesias  (2001) 
proposed  notions  of  time-varying  poles  that  rely  on  a  QR  factorization  of  the  state  transition 
matrix;  the  good  feature  of  this  notion  is  that  the  time-varying  modes  in  this  case  are  always 
bounded,  however  a  drawback  is  that  this  notion  does  not  necessarily  result  in  being  able  to 
decouple  the  modes  (i.e.  it  does  not  always  diagonalize  the  system  matrix,  even  though  it  may 
actually  be  possible  to  do  so).  Also,  these  poles  are  not  necessarily  related  to  the  LTI  poles  that 
occur  at  the  start/end  times  of  a  finite  duration  morphing. 

Using  Kamen’s  definition  of  time-varying  poles,  the  short  period  modes  of  the  morphing  MAV, 
for  two  different  morphing  trajectories  are  as  shown  in  Figures  46  and  47.  For  purposes  of 
comparison,  the  frozen  time  LTI  poles  are  also  shown  alongside.  The  feature  of  these  graphs  is 
that  when  the  morphing  occurs  from  +30  deg  to  -30  deg  over  a  span  of  2  seconds  (Figure  46),  the 
LTV  poles  remain  bounded  and  finite.  However,  when  the  direction  of  the  morphing  is  reversed 
(i.e.  the  same  MAV  morphs  from  -30  deg  to  +30  deg),  the  LTV  poles  become  infinite  (Figure 
47).  An  investigation  was  conducted  into  this  phenomenon  using  a  phase  plane  analysis  of  the 
nonlinear,  time-varying  Riccatti  equation.  Using  this  analysis,  it  is  possible  to  identify  precisely 
the  sets  of  initial  conditions  under  which  these  LTV  poles  (and  those  of  an  arbitrary  second  order 
time-varying  system)  diverge  to  oo.  Under  such  circumstances  an  optimization-based  approach 
in  order  to  determine  the  LTV  poles  using  only  output  measurements  was  explored.  The 
optimization  based  LTV  poles  are  found  to  remain  finite  under  all  circumstances  (even  when 
Kamen’s  LTV  poles  become  infinite);  additionally  they  compare  reasonably  well  with  Kamen’s 
LTV  poles  when  the  latter  is  also  finite. 
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Morphing  occurring  from  +30  deg  to  -30  deg 


Figure  46.  LTV  short  period  right  poles  (red)  and  frozen  time  LTI  poles  (blue)  shown 
while  morphing  occurs  from  +30  deg  to  -30  deg.  The  LTV  poles  remain  bounded 
throughout  the  duration  of  the  morphing. 


Figure  47.  LTV  short  period  right  poles  (red)  and  frozen  time  LTI  poles  (blue)  shown 
while  morphing  occurs  from  -30  deg  to  +30  deg.  The  real  part  of  the  LTV  poles  starts  to 
grow  unboundedly  after  t  =  0.7  sec. 
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COMPUTATION  OF  POLES  FOR  LTV  SYSTEMS 


The  fact  that  the  analysis  of  linear  time  varying  (LTV)  systems  is  substantially  different  from 
that  of  linear  time  invariant  (LTI)  systems  has  been  well  established  [6].  For  example, 
eigenvalues  of  the  system  matrix  (that  govern  the  stability  of  LTI  systems)  have  no  significance 
whatsoever  when  it  comes  to  the  stability  of  LTV  systems  [9].  Due  to  the  possibility  of 
nonexponential  responses  for  LTV  systems,  the  connections  between  asymptotic  stability  and 
bounded-input  bounded-output  stability  for  such  systems  are  not  as  straightforward  as  they  are 
for  LTI  systems  [5].  The  simple  notions  of  controllability/observability  in  LTI  systems  assume 
the  much  subtler  notions  of  total  controllability,  complete  controllability,  uniform 
controllability/observability  for  LTV  systems  [6].  In  short,  the  analysis  of  LTV  systems  requires 
considerably  enhanced  machinery  when  compared  to  that  used  in  the  analysis  of  LTI  systems. 


In  [4],  the  author  defines  a  notion  of  poles  and  zeros  of  LTV  systems.  For  an  LTI  system,  use  of 
the  Laplace  transform  converts  the  system  to  an  algebraic  representation  —  whose  numerator  and 
denominator  can  then  be  factorized  in  a  conventional  manner,  to  obtain  the  poles  and  zeros  of 
the  system.  For  an  LTV  system,  however,  use  of  the  Laplace  transform  does  not  (in  general) 
convert  the  system  to  an  algebraic  representation;  therefore  special  factorization  techniques  (as 
those  discussed  in  [8])  that  work  directly  on  the  ODEs  are  to  be  invoked.  Other  notions  of  poles 
and  zeroes  of  LTV  systems  have  also  been  discussed  in  [7],  [11],  [12],  In  [7],  the  authors 
compute  the  LTV  poles  subsequent  to  performing  a  QR  decomposition  of  the  state  transition 
matrix  of  the  system.  In  [II],  [12],  the  author  discusses  the  notion  of  Parallel  D  spectra  and 
Series  D  spectra. 


This  research  performed  a  detailed  phase-plane  analysis  on  the  differential  Riccatti  equation  [8] 
governing  the  evolution  of  the  poles  of  a  second  order  LTV  system.  The  objective  was  to 
identify  the  sets  of  initial  conditions  on  the  coefficients  of  the  LTV  system,  that  lead  to 
divergence  of  the  LTV  poles  to  negative  oo.  While  there  has  been  some  partial  work  in  this 
regard  in  [10],  that  phase  plane  analysis  only  considered  the  case  when  the  frozen  time  roots  of 
the  LTV  system  are  real.  Analysis  was  performed  under  this  research  project  under  the  different 
conditions  of  when  the  frozen  time  roots  of  the  LTV  system  are  wholly  real,  or  wholly  complex, 
or  transition  from  real  to  complex  at  some  time,  and  vice  versa.  Under  these  conditions  (that 
cause  the  LTV  poles  to  diverge),  we  then  explore  the  use  of  an  optimization-based  approach  that 
uses  output  measurements  to  compute  the  LTV  poles.  The  use  of  output  measurements  to 
compute  LTV  poles  would  be  necessary  in  practical  situations  (for  example,  while  performing 
flight  tests  to  detect  the  poles  of  a  time-varying  aircraft  with  uncertain  parameters).  The  ensuing 
optimization  based  LTV  poles  were  compared  with  the  LTV  poles  generated  from  the  Riccatti 
equation. 
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ANALYSIS 


Consider  a  second  order  LTV  system  of  the  form: 
d2x/dt2  +  ai  (t)dx/dt  +  ao(t)  x(t)  =  0  ( 1 ) 

which  can  be  written  using  operator  notation  D  =  d/dt,  as 
(D2+a,(t)D  +  ao(t))x(t)  =  0  (2) 

If  there  exist  functions  pi(t)  and  p2(t)  such  that  one  can  write 

(D2  +  *r(t)D  +  ao(t))x(t)  =  (D-p,(t)){(D-p2(t))x(t)}  (3) 

and  a  (non  commutative)  polynomial  multiplication  o  is  designed  such  that 

{(D-p,(t))o(D-p2(t))}x(t)  =  (D-p  i  (t))  ( (D-p2(t))x(t) }  (4) 

Then,  comparing  (3)  and  (4),  one  can  define 

D  o  p2(t)  =  p2(t)D  +  dp2/dt(t),  and  finally  arrive  at  an  equation  for  p2  as 

p22(t)  +  dp2/dt  +  a,  (t)p2(t)  +  ao(t)  =  0  (5) 

from  which  one  can  then  determine  pi  using 

Pi(t)  +  P2(t)  =  -ai(t)  (6) 

(pi(t),p2(t))  then  form  a  pole  set,  and  p2(t)  is  called  a  right  pole.  Note  that  this  is  an  ordered  pole 
set,  and  even  when  pi(t)  and  p2(t)  are  complex,  they  need  not,  in  general,  be  complex  conjugate 

[4]- 

The  mode  associated  with  a  right  pole  p2(t)  is  defined  as 
cpp2(t,0)  =  exp(j0l  p2(x)  dt)  (7) 

and  one  can  write 

x(t)  =  Ci  (pp2,(t,0)  +  C2  cpp22(t,0)  (8) 

where  p2i  and  p22  represent  two  right  poles  determined  by  solving  Eqn  (5)  from  two  different 
initial  conditions  (which  are  typically  the  two  time  invariant  poles  at  t  =  0,  also  referred  to  as  the 
"frozen-time"  roots  of  (1)  at  t=0);  and  Ci  and  C2  are  constants. 
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PHASE  PLANE  ANALYSIS  OF  POLES  OF  A  SECOND  ORDER  LTV  SYSTEM: 


In  this  section  a  phase  plane  analysis  of  (5)  is  performed,  when  used  to  compute  the  LTV  poles 
of  (1),  under  different  conditions  of  the  frozen  time  roots  of  (1).  Note  that  the  frozen-time  roots 
of  Eqn  (1)  are  given  by  {-at(t)  ±  sqrt(ai(t)2  -  4ao(t))}/2,  for  each  t.  These  are  also  the  same  as  the 
frozen-time  eigenvalues  of  the  matrix 

0  I 

-ao(t)  -ai(t) 

now  consider  different  cases  of  the  values  of  the  frozen-time  roots  : 

Case  I :  ai(t)2  -  4ao(t)  >  0  for  all  t,  i.e.  the  frozen-time  roots  are  real  for  all  time  t. 

In  this  case,  one  can  write  Eqn  (5)  as 

dp2/dt  =  -(p2  +  a,(t)/2)2  +  (a,(t)2/4  -  ao(t))  (9) 

Note  that  the  fact  that  the  frozen  time  roots  are  real  for  t=0,  and  that  these  same  roots  are  being 
used  as  initial  conditions  in  the  solution  of  Eqn  (9)  implies  that  p2i(t)  and  v,(t)  obtained  from  Eqn 
(9)  will  be  real  for  all  time.  One  can  therefore  perform  a  phase-plane  study  on  the  (p2,dp2/dt) 
plane.  On  this  plane,  at  each  time  t,  Eqn  (9)  represents  a  parabola  with  vertex  at  (-a i (t)/2 ,a i (t)2/4 
-  ao(t))  and  pointing  downward,  as  shown  in  Figure  48.  On  this  figure,  the  black  dots  represent 
the  frozen  time  roots,  while  the  arrows  indicate  the  direction  of  the  (p2,dp2/dt)  trajectory 
(determined  from  the  sign  of  dp2/dt  close  to  the  frozen  time  roots).  From  this  figure,  one  can  see 
that  once  the  left  branch  of  the  trajectory  enters  the  region  dp2/dt  <  0,  it  will  remain  in  the  region 
dp2/dt  <  0  for  all  future  time.  In  other  words,  once  the  LTV  pole  originating  at  the  left  (negative) 
fixed  point  begins  to  decrease  in  value  with  time,  it  will  continue  to  decrease  for  all  future  time. 
Since  the  y  co-ordinate  of  the  vertex  of  this  parabola  is  always  positive,  therefore  any  subsequent 
changes  in  ai(t)  and/or  ao(t)  (as  long  as  they  continue  to  satisfy  ai(t)2  -  4ao(t)  >  0)  will  not  alter 
this  behavior. 

Now  look  at  the  conditions  where  the  LTV  system  Eqn  (1)  will  have  trajectories  that  originate  at 
the  negative  initial  frozen  time  root,  and  continue  to  have  dp2/dt  <  0.  For  this,  look  at  the  second 
derivative  of  p2  and  evaluate  it  at  t=0.  Differentiating  Eqn  (5)  with  respect  to  time, 

d2p2/ dt2  =  -(2p2  +  ai(t))  dp2/dt  -dai/dt  p2  -  dao/dt  (10) 

The  above  equation  is  now  used  to  evaluate  d^/dt^O)  (keeping  in  mind  that  dp2/dt(0)=0),  to  get 
d2p2/dt2(0)  =  -da,/dt(0)  p2(0)  -  dao/dt(0)  (1 1) 
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Now  if  d2p2/dr(0)  <  0,  that  means  dp2/dt(At)  <  dp2/dt(0)s  (for  positive  At),  or  in  other  words, 
dp2/dt(At)  <  0  (since  dp2/dt(0)=0).  To  summarize  therefore,  if  the  parameters 
ai(0),dai/dt(0),ao(0),dao/dt(0)  in  Eqn  (1)  are  such  that  they  satisfy 

-dai/dt(0)  p2(0)  -  dao/dt(0)  <  0,  i.e. 

-dai/dt(0)  {(-a_{l}(0)-sqrt{a_{l}(0)2-4ao(0)})}/{2}  -  dao/dt(0)  <0  (12) 

then  the  LTV  pole  obtained  from  Eqn  (5)  and  originating  from  the  negative  initial  frozen  time 
root,  will  go  to  negative  oo.  On  the  other  hand,  if  the  parameters  satisfy 

-da,/dt(0)  {(-a_{l}(0)-\sqrt{a_{l}(0)2-4ao(0)})}/{2}  -\dao/dt(0)>0  (13) 

then  the  LTV  pole  obtained  from  Eqn  (5)  and  originating  from  the  negative  initial  frozen  time 
root,  will  remain  bounded.  In  the  case  when  the  inequality  sign  in  (12)  is  replaced  by  an  equality 
sign,  one  can  derive  corresponding  conditions  to  determine  whether  the  LTV  pole  is  bounded  or 
not,  by  looking  at  d'p2/dt3(0),  and  so  on. 

Now  each  of  the  above  two  behaviors  with  examples  will  be  illustrated.  Consider  a  system  with 
ai(t)  =  3-t;ao(t)  =  -0.1.  It  can  be  verified  that  these  system  parameters  satisfy  the  condition  (12). 
The  response  x(t)  of  Eqn  (1)  is  shown  in  Figure  49,  while  the  LTV  poles  computed  from  Eqn  (5) 
are  shown  in  Figure  50.  It  is  clearly  seen  that  within  the  time  frame  of  t=[0,5]  sec,  one  of  the 
LTV  poles  grows  unboundedly,  while  the  other  one  remains  bounded.  Figure  51  shows  why 
this  happens,  using  the  phase  plane  analysis  described  above.  This  figure  shows  the  phase 
portraits  at  several  time  instants  (as  aj(t)  varies),  along  with  the  corresponding  (p2,dp2/dt) 
trajectory  obtained  from  the  solution  of  Eqn  (5). 

On  the  other  hand,  consider  the  system  ai(t)  =  3+t;ao(t)  =  -0.1.  It  can  be  verified  that  these 
system  parameters  satisfy  the  condition  Eqn  (13).  The  LTV  poles  for  this  system  are  shown  in 
Figure  52.  It  is  seen  that  both  the  LTV  poles  p2|  and  p22  remain  bounded.  An  understanding  of 
why  this  is  so,  can  be  inferred  from  the  phase  portrait  of  Figure  (53). 

Case  II :  aj(t)2  -  4ao(t)  <  0  for  all  t,  i.e.  the  frozen-time  roots  are  complex  for  all  time  t. 

In  this  case,  p2i(0)  and  p22(0)  are  complex.  Separating  p2  into  its  real  and  imaginary  parts,  we 
can  write  p2  =  P2R  +  jp2i-  We  can  then  write  out  two  ODEs  from  Eqn  (5),  corresponding  to  the 
real  and  imaginary  parts  of  p2  as  : 

dp2R/dt  =  -p2R‘  +  P212  -  ai(t)  p2R  -ao(t)  (14) 

dp2i/dt  =  -2p2Rp2i  -  ai(t)  p2i  (15) 

Eqn  (14)  represents  a  hyperbolic  parabaloid  in  the  three-dimensional  (p2R,  P21,  dp2R/dt)  space. 
This  surface  is  as  shown  in  Figure  54.  In  general,  cross-sections  taken  parallel  to  two  of  the 
three  co-ordinate  planes  of  a  hyperbolic  parabaloid  are  parabolas,  while  that  taken  parallel  to  the 
third  co-ordinate  plane  is  a  hyperbola.  Similarly,  Eqn  (15)  represents  a  surface  that  is 
qualitatively  depicted  in  Figure  55.  The  solutions  of  Eqns  (14)  and  (15)  will  represent  curves 
that  lie  on  each  of  these  two  surfaces.  We  can  then  perform  an  analysis  by  qualitatively  plotting 
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the  trajectories  on  a  (p2R,p:i)s  plane.  Considering  for  a  moment  dp2R/dt  =  0  in  Eqn  (14),  then 
write  Eqn  (14)  as 

P21*  -  (P2R  +  ai(t))/2 )2  =  ao(t)  -  ai(t)2/4  (16) 

Eqn  (16)  then  represents  a  hyperbola  on  the  (p2R,  P21)  plane;  whose  two  vertices  at  each  time  t  are 
at  (-ai(t)/2,  ±  ao(t)  -  ai(t)2/4).  These  vertices  correspond  to  the  frozen  time  roots  of  Eqn  (1). 
Points  lying  "inside"  this  hyperbola  correspond  to  dp2R/dt  <  0,  while  points  lying  "outside"  the 
hyperbola  correspond  to  dp2R/dt  >  0.  Therefore  trajectories  that  lie  inside  the  hyperbola  have  a 
horizontal  component  to  the  right,  while  those  that  lie  outside  the  hyperbola  have  a  horizontal 
component  to  the  left.  One  can  perform  a  similar  analysis  on  Eqn  (6)  to  determine  that 
trajectories  that  satisfy  either  p2t  >  0;  p2R  <  -  ai(t)/2  or  P21  <  0;  p2R  >  -  ai(t)/2  have  a  vertical 
upward  component  (in  the  direction  of  increasing  P21),  while  trajectories  that  lie  in  either  P21  >  0; 
P2R  >  -  ai(t)/2  or  P21  <  0;  P2R  <  -  ai(t)/2  have  a  vertical  downward  component.  In  this  way,  after 
appropriately  partitioning  the  (p2R,  P21)  plane,  one  can  determine  the  direction  of  the  trajectories 
in  each  such  partition.  These  trajectories  are  qualitatively  depicted  in  Figure  56.  The  black  dots 
on  this  figure  represent  the  frozen  time  roots.  From  this  figure,  we  can  infer  that  as  long  as  the 
complex  part  of  the  LTV  poles  is  non-zero,  the  LTV  poles  remain  bounded.  However,  the 
instant  the  complex  part  of  the  LTV  poles  drops  to  zero,  it  then  remains  zero  for  all  future  time 
(this  can  also  be  seen  by  substituting  p2r=0  in  Eqn  (15)  and  noticing  that  it  leads  to  making 
dp2i/dt=0),  and  the  real  part  of  the  LTV  poles  subsequently  goes  to  negative  00. 

Now  the  case  when  the  frozen  time  roots  change  from  real  to  complex,  and  vice-versa  with  an 
example  wherein  ai(t)  =  3-t;  ao(t)  =  0.2  for  t=[0,5]  will  be  illustrated.  In  this  case,  the  frozen 
time  roots  are  initially  real  fort=[0,2.1]  sec,  then  become  complex  conjugate  for  t  =  [2. 1,3.9]  sec 
and  finally  become  real  again.  The  frozen  time  roots  are  depicted  in  Figure  57  and  the 
simulation  response  is  depicted  in  Figure  58. 

The  LTV  poles  for  this  system  are  shown  in  Figure  59.  In  this  case,  we  see  that  both  the  LTV 
poles  go  to  negative  00.  The  reason  for  this  can  be  inferred  from  the  analyses  discussed  above. 
P21  goes  to  negative  00  while  the  frozen  time  LTI  poles  are  still  real,  for  the  reasons  enumerated 
in  Case  I  above.  During  this  time,  P22  remains  bounded  (as  is  seen  in  Figure  60).  However,  the 
instant  the  frozen  time  roots  become  complex,  P22  finds  itself  on  the  real  axis  of  the  (p2R,p2i) 
plane,  and  from  the  analysis  of  Case  II  (and  Figure  57)  one  can  see  why  it  goes  to  negative  00. 
Subsequently,  even  though  both  the  frozen  time  roots  become  real  again,  P22  cannot  be  recovered 
again  to  have  a  finite  value. 

DETERMINATION  OF  LTV  POLES  FROM  OUTPUT  MEASUREMENTS  THROUGH 
AN  OPTIMIZATION  BASED  APPROACH: 

In  this  section,  the  determination  of  the  poles  of  a  LTV  system  from  output  measurements  will 
be  investigated.  In  particular,  those  cases  when  the  LTV  poles  are  determined  from  the  solution 
of  Eqn  (5)  go  to  00,  and  are  consequently  unable  to  provide  an  accurate  representation  of  the 
modes  of  the  system.  First  consider  the  case  when  only  a  single  measurement  source  (of  x(t))  is 
available,  from  which  the  LTV  poles  need  to  be  determined.  For  a  second  order  LTV  system  (of 
the  form  in  Eqn  (1)),  one  can  write  x(t)  in  terms  of  its  right  poles  p2i(t)  and  p2z(t)  as  : 
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x(t)  =  Ci  \exp(\int_0At  (p2i(t)  dt))  +  C2  \exp(\int_0At  (P22O)  dt))  (17) 

where  C\  and  C2  are  constants  that  depend  on  the  initial  conditions  and  the  values  of  the  LTV 
poles  at  t=0  (which  are  the  same  as  the  frozen-time  roots  of  Eqn  (1)  at  t=0). 

Discretize  time  with  equally  spaced  intervals  At  as  [ti,t2,...,tn]  and  then  write  out  the  above  as  an 
optimization  problem  involving  nonlinear  least  squares  minimization  as  follows  : 

Find  P2i(ti,t2,...,tn),  P22(ti,t2,....,tn)  such  that 

P2i(ti,t2,...,tn),  P22(ti,t2,.»,tn)  =  arg  min  I  K»l=n  [x(ti) 

-  Ci  eA{(p2i(0)  +  (p2i(ti)  + +  p2i(t,-i))/2  +  P2i(ti))  At} 

-  C2 eA{(p22(0)  +  (P22  (t,)  + +  p22(ti-,))/2  +  p22(t1))  At}]2  (18) 

where  the  trapezoidal  approximation  for  the  integral  in  Eqn  (17)  is  used. 

Figure  61  then  shows  the  LTV  poles  computed  by  solving  Eqn  (17)  for  the  case  when  ai(t)=3- 
t;ao  (t)=0.2.  It  is  seen  that  both  the  LTV  poles  remain  bounded.  Figure  62  shows  the  ensuing 
LTV  modes,  along  with  a  comparison  of  the  modes  obtained  from  the  solution  of  Eqn  (5).  It  can 
be  clearly  seen  that  these  modes  are  able  to  pick  up  the  unstable  segment  of  the  response  shown 
in  Figure  58,  while  the  same  could  not  be  detected  by  the  solution  of  Eqn  (5).  Furthermore,  even 
though  the  frozen  time  roots  in  this  case  change  from  real  to  complex  and  back  to  real,  Eqn  (1 8) 
generates  only  real  LTV  poles  throughout  owing  to  the  fact  that  Ci  and  C2  are  real  (which  in  turn 
is  because  the  initial  values  of  the  frozen  time  roots  are  real).  Figure  62  shows  a  comparison  of 
the  ensuing  modes  from  the  different  approaches.  While  it  is  seen  that  the  modes  corresponding 
to  the  solution  of  Eqn  (2)  cannot  demonstrate  the  unstable  segment  of  the  response  (which 
occurs  after  p2i  has  gone  to  negative  00),  the  modes  determined  from  the  optimization  approach 
are  able  to  demonstrate  the  same. 

Figure  63  then  shows  the  LTV  poles  computed  by  solving  Eqn  (18)  for  the  case  when  aj(t)=3- 
t;ao  (t)=-0.1.  It  is  seen  that  LTV  pole  P21  remain  bounded.  However,  a  comparison  of  P22 
obtained  from  the  solution  of  Eqn  (18)  with  that  obtained  from  the  solution  of  Eqn  (5)  shows 
some  mismatch  of  the  two  poles,  particularly  in  the  region  surrounding  t=3  sec  as  well  as  close 
to  t=5sec.  A  modified  optimization  function,  that  now  involves  the  use  of  two  measurement 
sources  (i.e.  x(t)  as  well  as  dx/dt(t),  is  used)  is  therefore  explored.  Since, 

dx/dt(t)  =  Ci  p2i(t)\exp(Jo‘  (p2i(t)  dt))  +  C2  P22  (t)\exp(V  (p22(t)  dt))  (19) 

therefore,  the  optimization  problem  takes  the  form  : 

Find  P2i(ti,t2,...,tn),  p22(ti,t2,...,tn)  such  that 

P21  (tl,t2,...,tn),  P22  (ti,t2,...,t„)= 

arg_{min}  \Sigma_{i=2}A{i=n}  [x(t_i)  -  Ci  exp{(  p2i(0)  +  (p2i(ti)  + . +  p2]  (tj.i))/2  +  p2i  (t,)) 

At} 
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-  C2  exp{(  p22(0)  +  (P22(ti)  + . +  P22  (ti-i))/2  +  P22  (ti))  At}]-  +  [dx/dt(t,)  -  C|  P2i(t,)exp{(  p2i(0) 

+  (P2i(ti)  + . •+  P21  (ti-0)/2  +  P21  (ti))  At} 

-  C2  P22(t|)exp{(  P22(0)  +  P22(tl)  + . +  P22  (ti-l))/2  +  P22  (ti))  At}]2  (20) 

With  this  objective  function,  it  is  now  seen  that  the  LTV  pole  P22  obtained  is  much  closer  to  that 
obtained  from  the  solution  of  Eqn  (5).  This  is  evident  from  Figure  63.  It  turns  out  that  this  new 
objective  function  does  not  change  P21  at  all,  from  the  previous  optimization  function.  This  P21 
however  is  bounded,  in  contrast  to  the  unbounded  P21  obtained  from  the  solution  of  Eqn  (5). 

CONCLUSIONS: 


A  detailed  phase-plane  analysis  on  the  differential  Riccatti  equation  governing  the  evolution  of 
the  poles  of  a  second  order  LTV  system  was  performed.  The  objective  was  to  identify  the  sets 
of  initial  conditions  on  the  coefficients  of  the  LTV  system,  that  lead  to  the  LTV  poles  diverging 
to  negative  00.  This  was  done  under  the  different  conditions  of  when  the  frozen  time  roots  of  the 
LTV  system  are  wholly  real,  or  wholly  complex,  or  change  from  real  to  complex  at  some  time, 
and  vice  versa.  Thus,  for  example,  when  the  frozen  time  roots  change  from  real  to  complex  and 
then  back  to  real,  both  the  LTV  poles  diverge  to  negative  00.  Under  these  conditions,  the  use  of 
an  optimization  based  approach  that  uses  output  measurements  to  compute  the  LTV  poles  was 
explored.  The  ensuing  optimization  based  LTV  poles  with  the  LTV  poles  generated  from  the 
Riccatti  equation  were  compared.  It  is  observed  that  both  the  optimization  based  LTV  poles 
remain  bounded;  with  one  of  these  poles  comparing  reasonably  well  with  the  stable  solution  of 
the  Riccatti  equation. 


0 

p2 


Fig  48.  Phase  plane  trajectory  of  (p2,dp2/dt)  when  frozen  time  roots  are  real 
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Figure  49.  Simulation  response  for  ai(t)=3-t;ao(t)=-0.1 


Figure  50.  LTV  poles  computed  for  ai(t)=3-t;ao(t)=-0.1  using  Eqn  (5) 
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Figure  5 1 .  Phase  plane  trajectory  of  (p2,dp2/dt)  at  several  time  instants  t=0,l ,2,3,4, 5  sec  for 

ai(t)=3-t;ao(t)=-0.1 


Figure  52.  LTV  poles  computed  for  ai(t)=3+t;ao(t)=-0.1  using  Eqn  (5) 
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Figure  53.  Phase  plane  trajectory  of  (p2,dp2/dt)  at  several  time  instants  t=0,l,2,3>4,5  sec  for 

ai(t)=3+t;ao(t)=-0.1 
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Fig  54.  (p2R,  p2i,  dp2R/dt)  phase  space 
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Figure  55.  (p2R,  P21,  dp2i/dt)  phase  space 


Figure  56.  Phase  plane  trajectory  of  (p2R,p2i)  when  frozen  time  roots  are  complex 
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Figure  57.  Simulation  response  for  ai(t)=3-t;ao(t)=0.2 


Figure  58.  LTV  poles  computed  for  ai(t)=3-t;ao(t)=0.2  using  Eqn  (5) 
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Figure  59.  Phase  plane  trajectory  of  (p2,dp2/dt)  at  t=0,l,2  sec  for  ai(t)=3-t;ao  (t)=0.2 


Figure  60.  LTV  poles  computed  for  ai(t)=3-t;ao  (t)=0.2  using  Eqns  (18),  (20) 
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Figure  61.  Simulation  response  for  ai(t)=3-t;ao(t)=0.2 


Figure  62.  Comparison  of  LTV  modes  computed  for  ai(t)=3-t;ao  (t)=0.2  using  system  parameters 

(5)  and  sensor  measurements  (Eqns  (18),  (20)) 
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Fig  63.  Comparison  of  LTV  poles  computed  for  ai(t)=3-t;ao(t)=-0.1  using  system  parameters 
(Eqn  (5))  and  sensor  measurements  (Eqns  (18),  (20)) 

MORPHING  WING  MAV  STUDY 

Vehicle 

A  vehicle  which  features  the  independent  multi-joint  capability  was  designed  by  retrofitting  an 
existing  aircraft.  The  basic  construction  uses  skeletal  members  of  a  prepregnated,  bi-directional 
carbon  fiber  weave  along  with  rip-stop  nylon.  The  fuselage  and  wings  are  entirely  constructed  of 
the  weave  while  the  tail  features  carbon  spars  covered  with  nylon.  The  resulting  structure  is 
durable  but  lightweight. 

The  wings  are  actually  separate  sections  that  are  connected  to  the  fuselage  and  each  other 
through  a  system  of  spars  and  joints.  These  joints,  as  shown  in  Figure  64,  are  representative  of  a 
shoulder  and  elbow  which  serve  to  vary  the  sweep  of  inboard  and  outboard.  The  range  of  motion 
admitted  by  these  joints  is  approximately  ±30  deg. 
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Figure  64.  Joints  on  Wing 

The  wing  surface  must  be  kept  continuous  for  any  configuration  of  sweeping  because  of 
aerodynamic  concerns.  This  vehicle  ensures  such  continuity  by  layering  feather-like  structures, 
as  shown  in  Figure  65,  within  the  joint.  These  structures  retract  onto  each  other  under  the  wing 
when  both  the  inboard  and  outboard  are  swept  back.  Conversely,  they  create  a  fan-like  cover 
across  the  ensuing  gap  when  the  inboard  is  swept  back  and  the  outboard  is  swept  forward.  The 
contraction  and  expansion  of  the  surface  area  created  by  these  structures  is  smoothly  maintained 
by  a  tract  system  implemented  on  the  outer  regions  of  each  member. 


Figure  65.  Feather-Like  Elements 

Spars,  formed  from  hollow  shafts  of  carbon  fiber,  are  placed  along  the  leading  edge  of  each 
wing.  These  spars  act  as  both  a  rigid  source  to  maintain  the  leading-edge  curvature  and  a 
connection  of  each  independent  wing  joint.  The  inboard  spar  is  translated  horizontally  by  a 
servo-driven  linear  actuator  located  inside  the  fuselage.  The  inboard  spar  is  then  connected  to  the 
inboard  wing  section  at  the  shoulder  joint  located  on  the  outside  of  the  fuselage.  The  inboard 
spar  then  connects  at  the  elbow  joint  to  outboard  spar  at  roughly  the  quarter  span  point.  The 
outboard  wing  region  is  activated  independently  of  the  inboard  region  by  means  of  a  servo 
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attached  at  the  elbow. 


Overall,  this  vehicle  is  able  to  achieve  a  wide  range  of  sweep  orientations.  Some  representative 
configurations  are  shown  in  Fig.  66  to  demonstrate  the  range. 


Figure  66.  Sweep  Configurations 


Modeling 

A  set  of  plant  models  were  generated  to  represent  the  flight  dynamics  at  each  configuration. 
Essentially,  the  forces  and  moments  affecting  the  vehicle  are  computed  using  an  assumption  of 
steady-state  conditions.  The  resulting  models  do  not  properly  relate  the  unsteady  aerodynamics 
but  will  include  the  dominant  steady-state  aerodynamics.  The  time-varying  inertias  are  then 
introduced  to  account  for  the  morphing. 

Also,  the  vehicle  is  constrained  in  this  example  to  simplify  the  configuration  space.  The  physical 
vehicle  has  4  degrees-of-freedom  in  that  each  inboard  and  outboard  and  sweep  independently  on 
the  left  and  right  sides;  however,  this  example  constrains  the  left  and  right  wings  to  symmetric 
sweep  with  the  outboard  having  twice  as  much  sweep  as  the  inboard.  This  constraint  limits  the 
system  to  a  single  degree-of-freedom  which  is  appropriate  for  the  longitudinal  nature  of  the 
altitude  change  associated  with  the  mission. 

The  flight  dynamics  for  each  configuration  are  trimmed  for  straight  and  level  flight.  Actually,  the 
thrust  is  held  constant  for  each  of  these  configurations  to  note  the  propulsion  is  not  affected  by 
the  morphing.  Such  an  approach  is  particularly  useful  for  relating  models  that  may  be  trimmed  at 
various  airspeeds.  In  this  case,  the  thrust  and  drag  were  held  constant  so  the  trim  routine  found 
the  correct  airspeed. 
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Simulations 


Simulations  were  conducted  for  fast  and  slow  varying  wing  sweeping  conditions  in  various 
mission  scenarios.  A  simulation  of  closed-loop  behavior  for  a  variable  wing-sweep  aircraft 
indicates  the  adverse  influence  that  such  transient  behavior  can  potentially  have  on  mission 
performance  if  not  accounted  for  in  the  flight  control  design  of  the  micro  air  vehicle.  In  this  case, 
the  coupling  during  the  transient  response  resulted  in  a  lateral  translation  that  caused  the  vehicle 
to  miss  its  target. 
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Task  5:  Computational  Aerodynamics  of  Flexible  and  Flapping  Wing  for  MAV 
MILESTONES/STATUS: 


First  9  Months/Status) 

•  Assess  computational  and  experimental  data  regarding  flexible,  morphing  wing 
aerodynamics  (lift  drag,  detailed  flow  field).  (Completed) 

•  Initiate  flexible.  Morphing  wing  simulations,  (completed  fro  flexible  and  flapping  wings) 

•  Document  and  evaluate  higher  order  numerical  and  modeling  treatments,  including 
numerical  dispersion/dissipation,  grid  movement,  and  transition-turbulence  model 
coupling  for  morphing  wings.  (Completed  for  flapping  wings) 


Under  the  leadership  of  Dr.  Wei  Shyy  from  the  University  of  Michigan  (UM)  and  Dr.  Jian  Tang ' 
post  doc  from  (UM)  significant  progress  was  made  in  computational  fluid-structure  interaction 
of  a  deformable  flapping  wing  for  micro  air  vehicle  applications.  Motivated  by  micro  air  vehicle 
applications,  a  fluid-structure  coupling  procedure  was  developed  between  a  Navier-Stokes  solver 
and  a  three- 
dimensional  FEM 
beam  solver  (  see 
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highlighting  some 
of  implications 
presented  in  this 
progress  report. 

The  fluid  model  _ 

includes  laminar, 
the  k-s 

turbulence  closure, 
and  a  filter-based 
k-e  closure.  The 
structural  model  is 
based  on  an  asymptotic 
approximation  to  the  equations 
of  elasticity.  Using  the  slenderness  as  the  small  parameter,  the  equations  are  decomposed  into 
two  independent  variational  problems,  corresponding  to  (i)  cross-sectional,  small-deformation 
and  (ii)  longitudinal,  large  deformation  analyses.  A  model  example  problem  corresponding  to  a 
NACA0012  wing  of  aspect  ratio  3  in  pure  heave  motion  is  presented  and  the  results  compared 


Fig  67:  Schematic  of  the  computational 
framework  for  flexible  wing  aerodynamics 
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against  available  experiment  data  (experimental  results  of  Heathcote  et  al1').  Quantitative 
comparisons  with  experiment  are  done  for  the  rigid  wing  and  the  implications  of  wing  flexibility 
on  aerodynamics  are  presented  in  a  qualitative  sense.  It  was  observed  that  phase  lag  of  the  wing 
tip  displacement  relative  to  the  flapping  motion  becomes  more  pronounced  as  the  fluid  density 
increases. 

SELECTED  RIGID  WING  RESULTS 

Figure  68  includes  the  thrust  coefficient  results  of  three  different  computational  cases  compared 
with  the  experiment  results.  It  may  be  observed  that  computational  results  have  a  decent 
agreement  with  those  of  experiment.  It  may  also  be  observed  that  the  inclusion  of  turbulent 
models  did  not  have  a  significant  effect  on  the  response.  The  experimental  data  shows 
asymmetric  thrust  in  the  two  strokes  of  each  cycle,  while  the  computational  results  show  no 
distinguished  difference  between  the  two  strokes.  Compared  to  the  experimental  data,  the 
computational  results  agree  with  the  experimental  data  very  well  at  the  small  force  peaks,  if  the 
asymmetric  difference  is  neglected. 


laminar 
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Fig  68:  Thrust  (upper)  and  Lift  (lower)  coefficients  of  experimental  and  computational 
results  of  rigid  wing,  k-e  model,  k-e  with  filter  model  and  laminar  flow  are  employed  in 
the  computation.  (Re=30,000,  hr=0.175,  k=1.82) 


SELECTED  FLEXIBLE  WING  RESULTS 

Fig.  69  includes  the  thrust  coefficient  response  of  the  flexible  wing  at  3  different  flow  densities 
(air  density,  lOxair  density,  and  416xair  density)  compared  to  the  computational  and 
experimentally  determined  responses  of  the  rigid  wing  at  the  same  Reynolds  number.  Once 
again,  the  effect  of  induced  phase  lag  is  seen  here.  Incremental  changes  to  the  thrust  coefficient 
are  clearly  seen  with  increasing  flow  density.  Also,  it  may  be  noted  that  spanwise  flexibility  has 
proven  to  be  beneficial  to  the  thrust  coefficient  just  as  seen  in  the  experiment. 
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Fig.  69  also  includes  the  experimental  response  of  the  flexible  wing  for  the  heavier  fluid  case. 
Similar  inferences  as  in  the  case  of  the  displacement  response  may  be  drawn  even  in  this  case. 


Fig  69:  Thrust  coefficient  response  of  the  flexible  wing  at  three  different  flow  dcnsitic 
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